{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampler statistics\n",
    "\n",
    "When checking for convergence or when debugging a badly behaving\n",
    "sampler, it is often helpful to take a closer look at what the\n",
    "sampler is doing. For this purpose some samplers export\n",
    "statistics for each generated sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "import pandas as pd\n",
    "import pymc3 as pm \n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a minimal example we sample from a standard normal distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pm.Model()\n",
    "with model:\n",
    "    mu1 = pm.Normal(\"mu1\", mu=0, sigma=1, shape=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [mu1]\n",
      "Sampling 2 chains: 100%|██████████| 6000/6000 [00:05<00:00, 1064.53draws/s]\n",
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/mkl_fft/_numpy_fft.py:1044: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  output = mkl_fft.rfftn_numpy(a, s, axes)\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step = pm.NUTS()\n",
    "    trace = pm.sample(2000, tune=1000, init=None, step=step, cores=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NUTS provides the following statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'depth',\n",
       " 'diverging',\n",
       " 'energy',\n",
       " 'energy_error',\n",
       " 'max_energy_error',\n",
       " 'mean_tree_accept',\n",
       " 'model_logp',\n",
       " 'step_size',\n",
       " 'step_size_bar',\n",
       " 'tree_size',\n",
       " 'tune'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.stat_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `mean_tree_accept`: The mean acceptance probability for the tree that generated this sample. The mean of these values across all samples but the burn-in should be approximately `target_accept` (the default for this is 0.8).\n",
    "- `diverging`: Whether the trajectory for this sample diverged. If there are many diverging samples, this usually indicates that a region of the posterior has high curvature. Reparametrization can often help, but you can also try to increase `target_accept` to something like 0.9 or 0.95.\n",
    "- `energy`: The energy at the point in phase-space where the sample was accepted. This can be used to identify posteriors with problematically long tails. See below for an example.\n",
    "- `energy_error`: The difference in energy between the start and the end of the trajectory. For a perfect integrator this would always be zero.\n",
    "- `max_energy_error`: The maximum difference in energy along the whole trajectory.\n",
    "- `depth`: The depth of the tree that was used to generate this sample\n",
    "- `tree_size`: The number of leafs of the sampling tree, when the sample was accepted. This is usually a bit less than $2 ^ \\text{depth}$. If the tree size is large, the sampler is using a lot of leapfrog steps to find the next sample. This can for example happen if there are strong correlations in the posterior, if the posterior has long tails, if there are regions of high curvature (\"funnels\"), or if the variance estimates in the mass matrix are inaccurate. Reparametrisation of the model or estimating the posterior variances from past samples might help.\n",
    "- `tune`: This is `True`, if step size adaptation was turned on when this sample was generated.\n",
    "- `step_size`: The step size used for this sample.\n",
    "- `step_size_bar`: The current best known step-size. After the tuning samples, the step size is set to this value. This should converge during tuning.\n",
    "- `model_logp`: The model log-likelihood for this sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the name of the statistic does not clash with the name of one of the variables, we can use indexing to get the values. The values for the chains will be concatenated.\n",
    "\n",
    "We can see that the step sizes converged after the 1000 tuning samples for both chains to about the same value. The first 2000 values are from chain 1, the second 2000 from chain 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c2124a8d0>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF79JREFUeJzt3X2QXfV93/H3BwlQwAEba+2hCFW4gRrhKgYrOB5ooFaMQXFM6CMMdkNLwtgOTEPoUDJmGEqcmZTU1PUMxgMdgo1qY8VpWiUoVRwH03Si2CwPEkg2toxt2OAxC8a1XdLw4G//OD/BZb3S3tW9u1e6vF8zd+45v3sevufs7v3sOb977klVIUnSQaMuQJK0fzAQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpWTrqAuZj+fLltWrVqlGXIUkHlHvvvffJqpqYa7oDKhBWrVrF5OTkqMuQpANKkm/1M52njCRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpKavQEhydpKHk+xKctUsr69McleS+5NsT7K+57U1SbYm2ZHkwSTLWvsX2jIfaI/XDW+zJEnzNed1CEmWADcC7wCmgHuSbKqqnT2TXQ1srKqbkqwGNgOrkiwFNgDvraptSV4LPNcz34VV5YUFOqA99t1n+Oy9U3g7Wi2ky9Ydz8FLFvakTj8Xpp0K7KqqRwCS3AGcC/QGQgFHtOEjgcfb8FnA9qraBlBVTw2jaGl/8qkvPcpNX/g6yagr0Tj7wD/6KQ5esrDr6CcQjgEe6xmfAt46Y5prgT9NchlwOPDzrf0EoJJsASaAO6rq+p75fi/JC8AfAB8q/8XSAehHPyqWHXwQX/mtc0ZdijSQfo4/Zvu/Z+Yb9wXAbVW1AlgP3J7kILrAOR24sD2fl2Rdm+fCqvoHwD9sj/fOuvLkkiSTSSanp6f7KFeStC/6CYQp4Nie8RW8dEpot4uBjQBVtRVYBixv895dVU9W1TN0fQuntOn+uj3/APgU3ampH1NVN1fV2qpaOzEx53czSZL2UT+BcA9wfJLjkhwCnA9smjHNo8A6gCQn0gXCNLAFWJPksNbBfAawM8nSJMvb9AcD7wIeGsYGSZL2zZx9CFX1fJJL6d7clwC3VtWOJNcBk1W1CbgCuCXJ5XSnky5q/QFPJ7mBLlQK2FxVdyY5HNjSwmAJ8GfALQuxgdJCs+NL46Kvr7+uqs10p3t6267pGd4JnLaHeTfQffS0t+3/Am+Zb7GSpIXjlcrSEGTWz15IBxYDQZIEGAiSpMZAkAbk9ZQaFwaCJAkwEKSh8HuMNA4MBEkSYCBIkhoDQRqQfcoaFwaCJAkwEKShsE9Z48BAkCQBBoIkqTEQpAHZp6xxYSBIkgADQRqKeKmyxoCBIEkCDARJUmMgSAPySmWNCwNBkgQYCNJQ2KWscWAgSJIAA0GS1BgI0oDKa5U1JvoKhCRnJ3k4ya4kV83y+sokdyW5P8n2JOt7XluTZGuSHUkeTLJsxrybkjw0+KZIkgaxdK4JkiwBbgTeAUwB9yTZVFU7eya7GthYVTclWQ1sBlYlWQpsAN5bVduSvBZ4rmfZ/xj44fA2RxoRe5U1Bvo5QjgV2FVVj1TVs8AdwLkzpingiDZ8JPB4Gz4L2F5V2wCq6qmqegEgyauA3wA+NNgmSJKGoZ9AOAZ4rGd8qrX1uhZ4T5IpuqODy1r7CUAl2ZLkviRX9szzW8CHgWf2pXBJ0nD1EwizHQzP7EW7ALitqlYA64HbkxxEd0rqdODC9nxeknVJ3gz8VFX94ZwrTy5JMplkcnp6uo9ypcXllcoaF/0EwhRwbM/4Cl46JbTbxcBGgKraCiwDlrd5766qJ6vqGbqjh1OAtwFvSfJN4H8DJyT5wmwrr6qbq2ptVa2dmJjod7skSfPUTyDcAxyf5LgkhwDnA5tmTPMosA4gyYl0gTANbAHWJDmsdTCfAeysqpuq6u9U1Sq6I4evVtWZw9ggaRTsU9Y4mPNTRlX1fJJL6d7clwC3VtWOJNcBk1W1CbgCuCXJ5XSnky6qqgKeTnIDXagUsLmq7lyojZEk7bs5AwGgqjbTne7pbbumZ3gncNoe5t1A99HTPS37m8Cb+qlDkrRwvFJZkgQYCJKkxkCQhsB7KmscGAiSJMBAkCQ1BoI0oPJSZY0JA0GSBBgI0lDYp6xxYCBIkgADQZLUGAjSgOxS1rgwECRJgIEgDYV9yhoHBoIkCTAQJEmNgSANyAuVNS4MBEkSYCBIQ+HXX2scGAiSJMBAkCQ1BoI0oPJaZY0JA0GSBBgI0lDYpaxxYCBIkoA+AyHJ2UkeTrIryVWzvL4yyV1J7k+yPcn6ntfWJNmaZEeSB5Msa+3/M8m21v7xJEuGt1mSpPmaMxDaG/WNwDnAauCCJKtnTHY1sLGqTgbOBz7W5l0KbADeV1UnAWcCz7V5/nlV/TTwJmAC+GcDb400Al6prHHRzxHCqcCuqnqkqp4F7gDOnTFNAUe04SOBx9vwWcD2qtoGUFVPVdULbfj7bZqlwCH4tfKSNFL9BMIxwGM941Otrde1wHuSTAGbgcta+wlAJdmS5L4kV/bOlGQL8ATwA+Cz8y9f2j94obLGQT+BMNuv+sz/5i8AbquqFcB64PYkB9H99386cGF7Pi/JuhcXUvVO4GjgUODts648uSTJZJLJ6enpPsqVJO2LfgJhCji2Z3wFL50S2u1iYCNAVW0FlgHL27x3V9WTVfUM3dHDKb0zVtX/Azbx46ehdr9+c1Wtraq1ExMTfZQrSdoX/QTCPcDxSY5Lcghdp/GmGdM8CqwDSHIiXSBMA1uANUkOax3MZwA7k7wqydFt+qV0RxVfGcYGSYvNzi+Ni6VzTVBVzye5lO7NfQlwa1XtSHIdMFlVm4ArgFuSXE7393FRVRXwdJIb6EKlgM1VdWeS1wObkhzalvnnwMcXYgMlSf2ZMxAAqmoz3eme3rZreoZ3AqftYd4NdB897W37DvAz8y1W2n/Zq6wDn1cqS5IAA0GS1BgI0oC8UlnjwkCQJAEGgjQUXqmscWAgSJIAA0GS1BgIkiTAQJCGwI8ZaTwYCNIQ2KescWAgSJIAA0GS1BgIkiTAQJAG5ldXaFwYCNIQeKWyxoGBIEkCDARJUmMgSJIAA0EamJ3KGhcGgjQE8VpljQEDQZIEGAiSpMZAkCQBBoI0sPLrrzUm+gqEJGcneTjJriRXzfL6yiR3Jbk/yfYk63teW5Nka5IdSR5MsizJYUnuTPKV1v47w9woabF5pbLGwZyBkGQJcCNwDrAauCDJ6hmTXQ1srKqTgfOBj7V5lwIbgPdV1UnAmcBzbZ7/WFVvBE4GTktyzuCbI0naV/0cIZwK7KqqR6rqWeAO4NwZ0xRwRBs+Eni8DZ8FbK+qbQBV9VRVvVBVz1TVXa3tWeA+YMVgmyJJGkQ/gXAM8FjP+FRr63Ut8J4kU8Bm4LLWfgJQSbYkuS/JlTMXnuTVwC8Cn59t5UkuSTKZZHJ6erqPciVJ+6KfQJjt7OjMXrQLgNuqagWwHrg9yUHAUuB04ML2fF6SdS8uuDul9Gngo1X1yGwrr6qbq2ptVa2dmJjoo1xpcXmlssZFP4EwBRzbM76Cl04J7XYxsBGgqrYCy4Dlbd67q+rJqnqG7ujhlJ75bga+VlUf2bfypf2DfcoaB/0Ewj3A8UmOS3IIXafxphnTPAqsA0hyIl0gTANbgDXtU0VLgTOAnW26D9H1N/z6MDZEkjSYOQOhqp4HLqV7c/8y3aeJdiS5Lsm722RXAL+aZBvdKaCLqvM0cANdqDwA3FdVdyZZAXyQ7lNL9yV5IMmvDH3rJEl9W9rPRFW1me50T2/bNT3DO4HT9jDvBrqPnva2TeFRtiTtV7xSWRqQfcoaFwaCJAkwEKShiN9doTFgIEiSAANBktQYCNKAvFJZ48JAkCQBBoIkqTEQJEmAgSBJagwEaUDeU1njwkCQJAEGgjQUXqiscWAgSJIAA0GS1BgI0qDsU9aYMBAkSYCBIA2FncoaBwaCJAkwECRJjYEgDcg+ZY0LA0GSBBgI0lAEe5V14OsrEJKcneThJLuSXDXL6yuT3JXk/iTbk6zveW1Nkq1JdiR5MMmy1v7bSR5L8sPhbY4kaV/NGQhJlgA3AucAq4ELkqyeMdnVwMaqOhk4H/hYm3cpsAF4X1WdBJwJPNfm+SPg1CFsgyRpCPo5QjgV2FVVj1TVs8AdwLkzpingiDZ8JPB4Gz4L2F5V2wCq6qmqeqEN/1VVfXvQDZBGrbypssZEP4FwDPBYz/hUa+t1LfCeJFPAZuCy1n4CUEm2JLkvyZUD1itJWiD9BMJsvWUz/yW6ALitqlYA64HbkxwELAVOBy5sz+clWTefApNckmQyyeT09PR8ZpUWjVcqaxz0EwhTwLE94yt46ZTQbhcDGwGqaiuwDFje5r27qp6sqmfojh5OmU+BVXVzVa2tqrUTExPzmVWSNA/9BMI9wPFJjktyCF2n8aYZ0zwKrANIciJdIEwDW4A1SQ5rHcxnADuHVbwkaXjmDISqeh64lO7N/ct0nybakeS6JO9uk10B/GqSbcCngYuq8zRwA12oPADcV1V3AiS5vvU5HJZkKsm1w944aTHYpaxxsbSfiapqM93pnt62a3qGdwKn7WHeDXQfPZ3ZfiVgJ7Mk7Se8UlkaAvuUNQ4MBEkSYCBIkhoDQRqQFyprXBgIkiTAQJCGIl6qrDFgIEiSAANBktQYCNKA7FPWuDAQJEmAgSANhV3KGgcGgiQJMBAkSY2BIA3IeyprXBgIkiTAQJCGw15ljQEDQZIEGAiSpMZAkAZkl7LGhYEgSQIMBGko7FPWODAQJEmAgSBJagwEaVD2KmtM9BUISc5O8nCSXUmumuX1lUnuSnJ/ku1J1ve8tibJ1iQ7kjyYZFlrf0sb35Xko/EehJI0UnMGQpIlwI3AOcBq4IIkq2dMdjWwsapOBs4HPtbmXQpsAN5XVScBZwLPtXluAi4Bjm+PswfdGGlU/H9G46CfI4RTgV1V9UhVPQvcAZw7Y5oCjmjDRwKPt+GzgO1VtQ2gqp6qqheSHA0cUVVbq/tmsE8CvzTgtkiSBtBPIBwDPNYzPtXael0LvCfJFLAZuKy1nwBUki1J7ktyZc8yp+ZYpiRpEfUTCLMdC8/sRrsAuK2qVgDrgduTHAQsBU4HLmzP5yVZ1+cyu5UnlySZTDI5PT3dR7nS4ip7lTUm+gmEKeDYnvEVvHRKaLeLgY0AVbUVWAYsb/PeXVVPVtUzdEcPp7T2FXMsk7a8m6tqbVWtnZiY6KNcSdK+6CcQ7gGOT3JckkPoOo03zZjmUWAdQJIT6QJhGtgCrElyWOtgPgPYWVXfBn6Q5Gfbp4v+JfA/hrJF0gjYpaxxsHSuCarq+SSX0r25LwFuraodSa4DJqtqE3AFcEuSy+lO/VzUOoufTnIDXagUsLmq7myLfj9wG/ATwJ+0hyRpROYMBICq2kx3uqe37Zqe4Z3AaXuYdwPdR09ntk8Cb5pPsZKkheOVytKAvKWyxoWBIEkCDARpKLxQWePAQJAkAQaCJKkxEKQB2amscWEgSJIAA0EainitssaAgSBJAgwESVJjIEiSAANBGpj3Q9C4MBCkIfBKZY0DA0GSBBgIkqTGQJAkAQaCNDC/ukLjwkCQJAEGgiSpMRAkSYCBIElqDARpQPYpa1wYCNIQxEuVNQYMBEkS0GcgJDk7ycNJdiW5apbXVya5K8n9SbYnWd/aVyX5myQPtMfHe+b5F23aHUmuH94mSZL2xdK5JkiyBLgReAcwBdyTZFNV7eyZ7GpgY1XdlGQ1sBlY1V77elW9ecYyXwv8LvCWqppO8okk66rq84NvkiRpX8wZCMCpwK6qegQgyR3AuUBvIBRwRBs+Enh8jmW+AfhqVU238T8D/gmwIIHwK5+4h2899cxCLFri8e/9DStfe/ioy5AG1k8gHAM81jM+Bbx1xjTXAn+a5DLgcODne147Lsn9wPeBq6vqL4BdwBuTrGrL+yXgkH2ovy8rjzqcQ5baXaKFcfzrX8WZf/91oy5DGlg/gTDbxydmftLuAuC2qvpwkrcBtyd5E/BtYGVVPZXkLcB/T3JSVT2d5P3AZ4AfAX9Jd9Tw4ytPLgEuAVi5cmVfGzXTNb+4ep/mk6RXkn7+bZ4Cju0ZX8GPnxK6GNgIUFVbgWXA8qr626p6qrXfC3wdOKGN/1FVvbWq3gY8DHxttpVX1c1Vtbaq1k5MTPS/ZZKkeeknEO4Bjk9yXJJDgPOBTTOmeRRYB5DkRLpAmE4y0TqlSfIG4Hhgd1/E69rza4APAP9l8M2RJO2rOU8ZVdXzSS4FtgBLgFurakeS64DJqtoEXAHckuRyutNJF1VVJfk54LokzwMvAO+rqu+2Rf/nJD/dhq+rqq8OedskSfOQOoC+zH3t2rU1OTk56jIk6YCS5N6qWjvXdH70RpIEGAiSpMZAkCQBBoIkqTmgOpWTTAPf2sfZlwNPDrGcYbGu+bGu+bGu+RnXuv5uVc15IdcBFQiDSDLZTy/7YrOu+bGu+bGu+Xml1+UpI0kSYCBIkppXUiDcPOoC9sC65se65se65ucVXdcrpg9BkrR3r6QjBEnSXox9IMx1P+hFWP83kzzY7ik92dqOSvK5JF9rz69p7Uny0Vbr9iSnDLmWW5M8keShnrZ515Lkl9v0X0vyywtU17VJ/rrnftzre177zVbXw0ne2dM+tJ91kmPbfcK/3O77/W9a+0j3117qGvX+WpbkS0m2tbr+fWs/LskX27Z/Jt03JpPk0Da+q72+aq56h1zXbUm+0bO/3tzaF+33vi1zSbp70f9xGx/p/qKqxvZB9+2sX6e7+c4hwDZg9SLX8E26e0P0tl0PXNWGrwL+QxteD/wJ3U2Jfhb44pBr+TngFOChfa0FOIruK8yPAl7Thl+zAHVdC/zbWaZd3X6OhwLHtZ/vkmH/rIGjgVPa8E8CX23rHun+2ktdo95fAV7Vhg8Gvtj2w0bg/Nb+ceD9bfgDwMfb8PnAZ/ZW7wLUdRvwT2eZftF+79tyfwP4FPDHbXyk+2vcjxBevB90VT0L7L4f9KidC3yiDX+C7haiu9s/WZ2/Al6d5OhhrbSq/hfw3RnN863lncDnquq7VfU08Dng7AWoa0/OBe6o7uZL36C7HeupDPlnXVXfrqr72vAPgC/T3U52pPtrL3XtyWLtr6qqH7bRg9ujgLcDn23tM/fX7v34WWBdkuyl3mHXtSeL9nufZAXwC7R7wbTtH+n+GvdAmO1+0Hv741kIRXe/6XvT3Q4U4PVV9W3o/sCB3TfkHUW9861lMWu8tB2237r71Mwo6mqH5yfT/Xe53+yvGXXBiPdXO/3xAPAE3Rvm14HvVdXzs6zjxfW31/8P8NrFqKuqdu+v32776z8lOXRmXTPWvxA/x48AV9LdRhi67R/p/hr3QOjnftAL7bSqOgU4B/i1dDcN2pP9od7d9lTLYtV4E/D3gDfT3Zv7w6OoK8mrgD8Afr2qvr+3SUdc18j3V1W9UFVvprvN7qnAiXtZx8jqSne/998E3gj8DN1poH+3mHUleRfwRHW3Fn6xeS/rWJS6xj0Q+rkf9IKqqsfb8xPAH9L9oXxn96mg9vxEm3wU9c63lkWpsaq+0/6QfwTcwkuHwYtWV5KD6d50/2tV/bfWPPL9NVtd+8P+2q2qvgd8ge4c/KuT7L4zY+86Xlx/e/1IutOGi1HX2e3UW1XV3wK/x+Lvr9OAdyf5Jt3purfTHTGMdn/ta+fDgfCgu0XoI3SdLbs7zk5axPUfDvxkz/Bf0p13/F1e3jF5fRv+BV7eofWlBahpFS/vvJ1XLXT/TX2DrmPtNW34qAWo6+ie4cvpzpMCnMTLO9EeoesgHerPum33J4GPzGgf6f7aS12j3l8TwKvb8E8AfwG8C/h9Xt5J+oE2/Gu8vJN0497qXYC6ju7Znx8BfmcUv/dt2WfyUqfyaPfXMDZof37QfWrgq3TnMz+4yOt+Q/thbQN27F4/3bm/zwNfa89H9fxy3thqfRBYO+R6Pk13OuE5uv8sLt6XWoB/Tdd5tQv4VwtU1+1tvduBTbz8De+Dra6HgXMW4mcNnE536L0deKA91o96f+2lrlHvrzXA/W39DwHX9PwNfKlt++8Dh7b2ZW18V3v9DXPVO+S6/rztr4eADbz0SaRF+73vWe6ZvBQII91fXqksSQLGvw9BktQnA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSAP8fZortzYRyFVgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(trace['step_size_bar'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `get_sampler_stats` method provides more control over which values should be returned, and it also works if the name of the statistic is the same as the name of one of the variables. We can use the `chains` option, to control values from which chain should be returned, or we can set `combine=False` to get the values for the individual chains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c212cc4a8>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXd8HNW5979nV5JluclFGFdkwBSbakQLJWCKQ0m4kEYSchNu8nKTm9yQvDc3F25CQt6ESwmBVEJIclOABG4ogUt3sI2pBskF9y5Ltrpk9bblvH9s0Wp3dndmd2Z3Vnq+n48+2p09c84zZ8785plnzjyjtNYIgiAIhYMn3wYIgiAI1hDhFgRBKDBEuAVBEAoMEW5BEIQCQ4RbEAShwBDhFgRBKDBEuAVBEAoMEW5BEIQCQ4RbEAShwChyotJZs2bpyspKJ6oWBEEYk9TU1LRprSvMlDUl3EqpbwBfBDSwGbhRaz2YrHxlZSXV1dVmqhYEQRAApdQBs2XThkqUUvOArwFVWuuTAC9wfebmCYIgCNlgNsZdBExUShUBZUCDcyYJgiAIqUgr3FrrQ8C9QB3QCHRprV9x2jBBEATBGDOhkunANcAiYC4wSSl1g0G5m5RS1Uqp6tbWVvstFQRBEABzoZJLgf1a61attQ94CvhAfCGt9UNa6yqtdVVFhakbo4IgCEIGmBHuOuAcpVSZUkoBlwDbnTVLEARBSIaZGPc64AlgPaGpgB7gIYftEgRBEJJgah631vp7wPcctkUQBEEwgTzyLgiCUGCIcAuCIBQYItyCIAgFhgi3IAhCgSHCLQiCUGCIcAuCIBQYItyCIAgFhgi3IAhCgSHCLQiCUGCIcAtjnoOH+znuOy+yu7kn36YIeeRLD9fwvWe25NsMWxDhFsY8L2xuZNgf5PH36vNtipBHXtraxB/fNv12MFcjwi0IglBgiHALgiAUGCLcwrhBqXxbIAj2IMItjHm0zrcFgmAvItzCuEGJyy2MEUS4hTGPONzCWEOEWxg3iL8tjBVEuAVBEAoMEW5BEIQCQ4RbGPPIrBJhrGFKuJVS5UqpJ5RSO5RS25VS5zptmCDYjgS5hTFCkclyPwVe0lp/TClVApQ5aJMg2IqWeSXCGCOtcCulpgIXAp8H0FoPA8POmiUI9qPE5RbGCGZCJUcDrcDvlVIblFK/VUpNctiulKy4fy33r9yVsHzIH+DE217i2U0NLP/xGn6xanf0t3PvfJXfvbE/l2aOoratj8pbnueE215MWe7x9+o45faXCQaz8xI//Zt3+M7fNhv+9uaeNipveZ7Ofnedfzv6hqm85Xne2dee8NuWQ11U3vJ8wt/G+s6EsoGg5uTbX+aJmoOA/THumx/bwMm3vxztwy2HujjvrlVsMrDFCo+uO8AZP1gJwKX3vcbPXt2dZo3M+fe/buLG37+bssxLW5o47jsvMjAcyLid+o5+Km95npv+VJ227IOv7eWCe1YB8Nbe0Bjt6AuN0b2tvdF9vqHusCUbWroHqbzl+YTly36wkkfeSZ4tsLVniMpbnufaB95M28Znf7eO/3za+HhzAjPCXQQsA36ltT4d6ANuiS+klLpJKVWtlKpubW212czR7Gzu4acGg7q9d5gBX4A7X9jOvtY+7n1lRNwbuwb5wXPbHLUrFS9uaQJg0BdMWe47f9tC96AfXzB1uXS8tbedR96pM/ztV2v2ArD5UFdWbdjN+gOhA/I3a/cl/Pbnd4235WGDNJ39w356Bv0JuZftenDymY0N9Az6gVAfrt7RwqHOAVZua86q3m8/vYX2vmG01uxp6eU+A+fELv5ac5DVO1Mfp/e8vINhf5BDnQMZt/Pq9lCfvGKib+56cQf1HaG2fv1aaAxsOhg6GT67sSFa7o9v1Vqy4W0DRwBCjsJ3/pY8P3fEgdhQl/6E/PruNv68zniMOoEZ4T4IHNRarwt/f4KQkI9Ca/2Q1rpKa11VUVFhp42CkDTI4TH4wRNW6CwvWkyhULZH0HNhtxkiVypuyBSgk3zOuD4Tl2FBF09HSivcWusmoF4pdXx40SVA/lzXMY7EYY1JJh4egx8ii3J1U9JugQu4RLkj4uaKERkjolb11Ki8izXZFGZnlfwr8Gh4Rsk+4EbnTBrfyAwIY5Kd0IzEMlI2Jx53TPt2CZxbPL2IFW5IzpVNj2R6TLllPxhhSri11huBKodtGdO4YOwXNEYhETAWlciBGn857MQusLNOjwqdbNwiGBEzkvV9LontEltCJSbKZHmbyVHkyUmX4pJj1zUk8/qMRCXSd7mKONh1lRTZRreESoLRUIkLlDsGM/Hp0eUzq8MtJ1AjRLhdhovHSl5JdsWS6komcuBZPdAtEdt+lpdVkbVdotsuuzkZE+O2um6G/enmY1GEO0e4YOwXNMm8PqObk5HjLf7Ac0KA7PRGo7Nh3KLcLmLUvrShe0yFSlys3CLcQkFgZVZJrrHt+A5vilsEwy12QPx0QIuhkgzbdPP5U4TbpbjomHEFyW9OJi5zNDRi0H509kW2dYX/B1yy811iBhB3c9LydMDEFczU4aYTVzwi3C5FpgWOJvnNyeShkuj3wghxx4RKsqvHLiJj0A0CNirGbUuoJH0luXQArCLC7TLcO1TySzJNtKKVbpsdEY9yXagk9N8V5mRhg4RKhIyx6o254mBxEUk9boMYSnzfOdmVSinbdlY0VOISxYhsljusGcHy1WiGT0665QRqxLgQbjdf8iSj8Cx2lqTTAY0WJuk8R2aVjHpyMsvpgOHK3DNc3RQqifmcI3Nccv40RDkhalVVVbq6On0ax3jiUy8unFHG3R89hU/95p3osns+egrfevJ9AGrvuooP/WQtZy+awSfOXMBVP3sjoc7rTp9H54CPVTtaEn67+PgKVu9sZV75RCpnlfH+wa5o1rcIV58yh+feb2TBjInUdwxw+ZLZCZnOjqmYREmRl3+77Di++Kdqli0sZ+ncaTwckzLyw6fO5X83jWQ4e/CGZfQOBbjrxR28cPP5nHXHq4Z9cvzsKTR0DrD5+ytY/uM17Gvt469fOpczK2dwzH++EPXOPnXWAm67egnL732Nr12yeFSKycqZZdS29/P2rcuZM21itJ8f/sJZlBZ7+fiDb3Pd6fN4asOhUW0/+eUP8NFfvZVg04qls3l5azPP/ev53PH8dmZOLqGxa5CaA4eZMamEy06czePV9WnribWvd8hP94Cf4UBigDfd+hGuOOlI7v/kaZxw20spyz3wmWX8y6PrAbhu2Tzu+8RpfOnhGlZub07p7X7/I0u5dMlszrtrVdIytXddZZhC1AqxY+ysRTOomDKBX356GRffu4b9bX0p1139zYu44qdr+dDSI3lzbzu//uwZLFs4nbr2fi780eqU686aPIFnv3oec8sn4g8EOfbbiSmIy0q89BukeV22sJz1cVn0fnr9adz82Mbo9zf+42LOv3s193/yVL7x+CZg5GnReC48roK1u5JnL5w1eQLXn7mAP75Vy+bvrwDgK4+up2vAx77WXhq6BlNuq1185eJjePy9g7T1DgGh/Z8pSqkarbWpJ9TN5irJC3Ud/fx67d5Ry+5+aceo7zuaetjR1EOR1/jiIV6MYomktTzUOZA0deVz7zcCRNNNGqWn3NsaOph+uWYPAOvrOhMGcaxoA/xi9R66Bny09Q6x+WDy9Ko7m3uin/eF23nknQOcWTljlMj85d16vrp8MU3dg/zo5dF9VNveD4TSVF57+vxRv/3hzVrAuJ/uW7nT0KaXt4b64JmNhxJSZnb0DY8SbYB7XzauJ96+ZBjlXjfixS1N3P4RX9pysXmun1p/iPs+cRpv7GlLG6K444XtTCzxmrIlG2LH2Lv7OwD45adJK9oQSqM66Avyt3Aa1DU7Wli2cDqrdyY6LvG09Q6xr7WPueUT6UuSg9tItIGE8Q6MEm2Av4e3KyLakNyrTTeXva13iF+s3jNq2fObG1Ou4wS/XL03fSEHcH2oJP5gShbrdMPlnFWcmj3gDxj3hd1dZHYOdbbT26yEODJtykwTbok9pyL+OIgcL2aPj2hfO7CpRvcjhMwoOOFOtu8LULdH0mZmMJ6N56amjkna3Udms8Zl+ySglYdszNy4Mk5MNTaIjzRlmpvciemodu/H8YzrhTvee0y2793gcecyJ0KqHMPJDlK7e8jscZhLj9uMQGXq9xXCTe5Ejzv0347ETNliSbjd39V5xf3CHRdPSLbz3SDcmZKJ6am2N7nHbeClZyHnZq98s40wWDvgC3cc2EH81U1kH5kOlYT/O9GLViIl43w3psX1wh0frk12ELthR1ueq22lbNwGGolhJvNus+k3s4KarZjm4oA3Y6OZ0FC+TxzxVzdummLoVKgk332eD9wv3Cbv4LnhvlHGopFB3UYDO9kLBFK1k82YN3sYZn9Tz95LbCfzUuVbQ+K7euRpTGv1OCGGTt1kznef5wPXC3d8jNuT1OL8771MLTBzkMSXSBXjTlqd0ToprE5nlumbk1mHSsyXzfdNrXyPwsRQSdjjNmtZ9H2d9mPN4zZPvvs8H7heuBOmAybxvlyRmMfBG0CJoRIjjzv5b6HfjWaimLchHrPHYbazSnIxHdDcVY+JE6zbQiXh/3a8YDdbLF3pZHFsjAcKTriTeV+FeHPSincYr32ptDDprJICncdtxVMzMw6MU8FasShV+/bUk3n7xh632ZPnyFVbfqcDWiHffZ4PXC/c/njhTqLcbth5TpoQL/Kp5nEnr8PcsgjpjjOvyRhG1vO4LcRKzLSUaU4RUzcn83zhnhAq8Vibx+1kYilLV05Wbk6Ow2CJaeFWSnmVUhuUUs85aVA8iaESY9xwuZTLm5OGs0os1pFsmd1kezVk6Qo7z8Mg3+0n3JyMLjfpcTuYWMrsPRGQm5PpsOJx3wxsd8qQZCSGSpJMB8yFMWmw/EqlLIzO9K0eBjVlbIPZS99czuPOdHvGiteWLLRoduvSPcSVK8bG3nAOU8KtlJoPXAX81llzEimkB3CcfHIyfvuM10kXKrH35qTZCEa20wGtzCox05STIp3vcZgQ4/ZE5nGbsyuyviPTAS2UtdJ+vvs8H5hK66qUegK4E5gCfFNrfXWq8naldRXGBnOmldKYozSbY5Xvf2Qp33t2a07amjOtlL/8n3O46N41ttb7y08v4yt/Xm9rnRFKijwM+/M/tSxXaV3TetxKqauBFq11TZpyNymlqpVS1a2tyfPoCuOPcegQ2c49cemMnaSxa5ABn3H61mxwMjmgG0Q7l5gJlZwHfEQpVQs8BixXSj0SX0hr/ZDWukprXVVRUWGzmUIhM1bix+MJJ8IPMgrsI61wa61v1VrP11pXAtcDq7TWNzhumTBmyPeNLsE6TjzQJlde9uH6edxC4eOGqZqCNeInBdjBeLyJ6BSWXl2mtV4DrHHEEkEQXIMTIivCbR/icQuOI6GS7Ml1Fzqxz0S37UOEW3AcCZUUHk68X1M8bvsQ4RYcRw7X7Mn1a3adEW7bqxy3iHALjiOOVvbkugvF43Y3ItyC48gBW3hkm4rXEBkGtiHCLTiO6HbhkW0qXiPkQSz7EOEWhAJAYtxCLCLcguNIqCR7cj8dUGLcbkaEW3AcOWALj/g3T9nBePC4fYHcJLsS4RYcZ9A3vjK3OUH/sP3Z+lLx1T9vsL3O2/62xfY63cbib7+Yk3ZEuAVBEAoMEW5BEIQCQ4RbEAShwBDhFgRBKDBEuAXbsPQydkEQMkaEW7ANryi3IOQEEW7BNjxOvg1WEIQoItyCbYjHLQi5QYRbsI0i8bgFISeIcAu2IaESQcgNItyCbYhuC0JuSCvcSqkFSqnVSqntSqmtSqmbc2GYUHh4RbkFIScUmSjjB/5Na71eKTUFqFFKrdRab3PYNqHAEOEWhNyQ1uPWWjdqrdeHP/cA24F5ThsmFB4q5+n+BWF8YinGrZSqBE4H1jlhjFDYLJk7Nd8mCMK4wEyoBACl1GTgSeDrWutug99vAm4CWLhwYUbGPP0vH+Abj2+ko28YgJPmTWNCkYebLjyGB9bswetRdA/4WF/XyZI5U6nr6GfGpBKKPIrJpUVUHTWDA+19KKXY09LDhCIvR0ydwMBwgEuXzGZ6WTHv7Ougb8jPxBIv9R390TzHFyyeRX3HAG/va+e0BeWctWgG6/aHytYcOAzAiqWzmTaxmEFfkJ5BH16Ph89/oJKnNxyivW+IshIv5x4zi9+9vo/a9n4uXzKbAV+APS29nL1oBpWzJvHB4yr4/O/f45ITjsAX1MyfPpFdTT3saunhipPm0NQ1SM+gj7MWzeRQZz/9QwG6B310D/o5fvYUag4c5uITKphaWszull4mTyiKbvP86RNp6hrEFwjS0jPEiXOmcuqCclZtb0YDR8+aTFBr2vuG+d9NDXz36iXsbe3lufcbqZgygQsXV3DuMTPZ1tDN4+/VsWBGGXUd/cyeWsrc8lImTyhia0M3XQM+zjtmFtubQp9POHIK5x49k2tOm8etT21mwYyJbG3oZtnC6by7vwON5qLjj2D21FJe2NzInGmlFHkU9YcH2NfaS+eAj//40Ak8uu4AtW39LD/hCDbWd3L+4lkM+gLUtffzzx88hiOmTODB1/ZyoL2fylllfHTZfH6+ag/+YJDlJ8zmgsWz2HKoi/1tfRzu93HS3KnMn17GqzuamT+9jKc3HOTkedN4Y3cb86aX8YNrlvLoujqWzg2NpUtOnM3/vFePUjCx2AvAiXOm8uymBr52ybE8s7EBj1K8vLWJKaVFeJXiXy4+lq0NXSydO43n3m8A4Mippby6vYXFsyez/IQjuHzpkTz+Xj1v7G6jc2CYMytnsK+1j22N3VTOLOPk+eUcaO/jshNn09o7RF1HP5895yjqO/p5Yv1B+ocDlHg9HHvEZMpKvEyeUMzMySW09gzR2DXAwcMDXHz8EZy6oJy1u1oZ8AU4bUE5h/uG+fXafUyfVMyQL8iVJ89hYomXYFDjC2iCWrO3tZf50yeyt6WPT529ILqNsyaX0D8c4PjZUwhozbFHTCYQ1Kze0cKmg12cNG8a/UN+JpcW0dw9yLzyMqoqp7N2VyulxV6eez9Uz8xJJcwtn0h5WQmbD3XyqbMWsrOph/1tfWgNXQM+Zk8r5cOnzOGp9YeYPKGIjv5hjpxaSnP3IB8+dS51Hf109vs4e9EM5k+fyANr9jLsD3L+4lmsWHokG+oO88LmRna3hLbFH9CUFHk4fWE5JV4PXQM+egb9bKjvZPbUUo6pmMRxs6fw9+3NnLagnDf2tHHu0TMpKfJQ7PXw5p42/uG0eRR5FfUdAwSCQcomFHHo8ADLTziCldua2dfWR117H33DAc6snI5HKRbPnsxnz6nMSPusorSJt5MopYqB54CXtdb3pStfVVWlq6urbTBPEARhfKCUqtFaV5kpa2ZWiQJ+B2w3I9qCIAiCs5iJcZ8HfBZYrpTaGP670mG7BEEQhCSYCpVYrlSpVuBAhqvPAtpsNMcuxC5riF3WELusMRbtOkprXWGmoCPCnQ1KqWqzcZ5cInZZQ+yyhthljfFulzzyLgiCUGCIcAuCIBQYbhTuh/JtQBLELmuIXdYQu6wxru1yXYxbEARBSI0bPW5BEAQhBSLcgiAIBYYItyAIQoEhwi0IglBgiHALgiAUGCLcgiAIBYYItyAIQoEhwi0IglBgiHALgiAUGCLcgiAIBYYItyAIQoEhwi0IglBgiHALgiAUGCLcgiAIBUaRE5XOmjVLV1ZWOlG1IAjCmKSmpqbN7DsnTQm3UuobwBcBDWwGbtRaDyYrX1lZSXV1tZmqBUEQBEApZfoF62lDJUqpecDXgCqt9UmAF7g+c/MEQRCEbDAb4y4CJiqlioAyoME5k9zB3tZefIFgvs0wpL6jn/5hf1Z1DPkD7G/rS1tud3MPwaCmd8jPwcP9WbUZy6AvQK2J9p3A7m1xE2/tbaOtdygvbTd1DdLV7xu1rLl7kM7+YQC01uxq7jFcN9Vv+aar30dTV9IAQ15IK9xa60PAvUAd0Ah0aa1fiS+nlLpJKVWtlKpubW2139Ic0tg1wCU/fo0fPrct36YYcsE9q7nht+uyquOWJzdz8b1r6B70JS2z+WAXl92/ll+v3cd1D7zJ+XevzqrNWG5+bAMX3buGIX/AtjrN8vEH37Z1W9zC67tb+fRv1lH1w7/npf1z7nyV8+5eNWrZ2f/1Kst+sBKA/36zlsvvX0vNgcMJ6/4+/Ft1bUdObLXCBfes4pw7X823GaMwEyqZDlwDLALmApOUUjfEl9NaP6S1rtJaV1VUmIqvu5aOvpCHsG6/+wZRhPV1nVmt/8aeNgAGh5MLZ33YK91U38mu5t6s2otnzc7QyT2Yh4ua7Y3duW80B+xtsXcfZULvUOKVYDD8WttN9aExW9+ReLWz+VAXAAfa3Xcl1D2Y3dWtE5gJlVwK7Ndat2qtfcBTwAecNUtwE0rl2wLBDKpAdlSBmOlqzAh3HXCOUqpMhUbGJcB2Z80S3IDWDtYd/e9gI4IwRjET414HPAGsJzQV0AM85LBdgotwxEMK63VQdNs23O7JptrVLjfddZiax621/h7wPYdtEcYREU9bO+nWjzPcLn6yr+1DHnlPQaHEDJ3CyTCGFo973JLquJLhYA4RbiEtykFfLp9e2JjzAMe5ozGeEOEWkpKLm5P59LgD4u7nFOlt+xDhNmCsOWJZ44AjF/F28+n1jjXdLhR/29DOQjHeJYhwC0lxUtfc4HEH5QwtFCgi3AZIqDBExBt2ZDagHt1GPhhroRLXj9ux1d15RYTbAHHEcod43Pbh5E1kO4jMUkp1ghlzN4wdQoRbyCv5fHIyH3lSnMT1HncYoxOM2086bsP1wr2h7jBdA6EMdkP+ABvrO9nXai6ZTs+gL+vUoYGgZltDN3taeqg5kDrpVFvvEI1dAynLbG3oIhjjZmqt2RJOsJOMA+190T7Y2jBSdndzD1prNtV3sqMplDipsWsgmtbTVN0d/TR0DkT7qb6jnzd2t7GtoTvqkaaad9vUNUhrzxAdfcMc6hzZ9r2tvXT0DbOxvpMD7SP7YMuhrlFe1ZqdrfSkyFBY39EfTQtqhZaeQZq7U6fijPe4g0E9qn9j2dnUw7B/tNLvbu5h0BeguXuQlp5QW7F9PuwP8s6+9lH7O8K2hu5oqKa2rY/9bX0ZpQ49eHikf2L3UiTtb21b36j+7R3y8/LWprTjNB3dg75R+zWW1p4hW9KgRvo3nv5hP3tNaED8WAPY0zJS55ZDXbT1DjHsD7KzyVpK2c7+YcNkWbnCkVeX2UXXgI9rH3iLFUtn8+vPVvGfT23hyfUHAai966q063/8wbfZ0dRjqqwRCvj5qt385O+7o8tWfuNCFs+eYlg+kk4zWXsb6g5z7QNv8e8rjucrFx8LwCPr6rjtb1v40z+dxYXHGWdV/OCP1nDUzDJ+/PFT+diDb0eXX3b/Wm6+ZDE/fTVk32v/fhEf/NGaqA0Pv3OA7z6zlUe+cDbnL55lWPfHY+qrvesqLrhnJN3pUTPLDNeJJZLustir8AU0tXddRTCoueTHr40qV3vXVaze0cKNf3iPO687Obr81qc288g7B3j+axcY1n/BPauZOamEmtsuS2tLLGfd8Wq03WQE4g7q37y+jztf3MGTXz6XM46aEV3e3D3Iip+s5VNnLeDO604B4HDfMJfdv5Z/OG0uf9vYEG3rrzUH+dYT7/Pbf6yitr2PHz6/nd99ropLTpwdrW9HUzdX/ux1vnrxsXxzxfFcdO+a6G9Wx+r5d69mSmkRm29fMWr5F/5QzV9uOoeL7l3DkjlTeeHmUP/e8Nt1bKzvZNGsSaz+5kWW2orlugfeYk9Lr6G9Z95hfByYiYJEinQN+Ljs/rV8+NS5/PxTp48q88U/VvPW3vaUfbVqRzP/9Idq7v7oyXzyzIVA6KR16X1rufLkI3ngM2dw9c/f4PSF5SydO5VH3qnjrVuWM7d8YnojgeU/fo2OvuGMtSVbXO1xRzycV7e3APCexVy9OyyeRY3YfHC0B9bSk3mS+obOkBcS69XtCKcYPZDm7H2gvT+aZjWW6pirgPa+0Z7p9sbQ9tcm8YzSEUmxaeYi1heIuYpIUmZf2KuP9262NqROsxq/XXYRLySR1KIHD4/2RiNXO9W1I3mk+8Ie7btxqX93hbdtf1sfu8OpcFvjxkzEG30/zdWQWXrCaUdjL4ze3tce/bwtJo3txnBqVTMv0UjFnixSyJoJ6USuGN4zSK381t72hGXx7GsNbV+sBkQ87Xf2dUQ98Q11ndH92tmf/Movng6HxqRZXC3c8fHPfMRD4weZHTe07Lz/EltXfL0eFVmeXYOFEju1SrLxZCbVQaSMmZ6NLxP5Pka71RYiQzbTsZdqHyqMj8FCGueuFm43ko0GZjsw0rc9uoAn3KAbZ7254iDRKb+mxIz5SbfRof3h9ht8qcZvfF/Z1UXJHBsXHhKWcLdwu6J3R4+oQvS4cz3tzYyH7waJyaZXVPRqxkQ7Scq44uSVQ6LTAS2sk6psqnGWar2Ek4QrdMYa7hZuF5LVwZ5t22kaj/9Zudrjzr9qZXX1RCRUkkI8kmxiJgJm0qiCIOWuj3swy85xEruvYkXfzPxyt+Fq4U51XFmJ2xpNxzKDUkZnZxeqYJh407xhl9uszcnKWR3P6VpTyh0akzTGbaWODIbDSPzW3l6Irc2NIpQyVGKhrKn6Um6/Mtzzbg81xeJq4Y4nVVggFXaGCux4aCPTm6zp1ooXXquhklx65h4XKEtCt1jYfivmJ9xkd++5P0fkbt/rJIHtQt8HrhbuVJ1rpd/j5+taIcETyLgm572geNus3pxMlrsj1jM0472b6u486nY0Pp3mdzOk3lTjipyaVRK7n/J/WswOM7NKUvV9uhj36LCJJdNcgauFO4JRv1oJWdi5Y/KZ38JqmGYkxm3W404fOshm82Ptd4OwxPenlSuhbOwfid9mUUka3HAPIRvsijvrZJ9lOqBzpLrEtORxZxEDsPcOdCTmnNnahiewFJ7DyDxuc/WbEXhz85ZNzCrJp8cd/p90tocFWTYTZ01exsEYt60120Mmwz7Vvkiv4NTkAAAVFUlEQVQ5q0QlHmtRLz4DO9yGu4U7RUdbET+roZLY4okDJ/uTgFM+e2KMO+xxmzxxJT3BxXRBpqESrXXcpXz+Dh+jg9oq2uBTQjtJShTglbmtpA5/hHrHyuPxVttQqvDndJsSbqVUuVLqCaXUDqXUdqXUuU4bFku2HattzAKX16l1hoKYvPjIzUlz1Tu5bVrHhUrc4HHbcMhmN6sk6+ZHEVufGy/7M3kAJ/tQSeK0v0S7HJqe6SBmk0z9FHhJa/0xpVQJkD77kA2knA5o4aDL6uakA5P1bX0AJ8lnyCDGHafcHpUo5pmanmBbhvXYgZWHZ5IRfRGElcLxdmTevHF9scLtShmyLpCpH8DJbD2FMtQPN57skpFWuJVSU4ELgc8DaK2HAccyrDR2DaA1zC2fyJu724DQJfyu5p5RaUP/5716Pn32UXg9il3NPTR3D3LBYuPseoO+AK09QwwMB5g+qZhVO1pYsfRIAkHNL1bv4QvnL2LW5AnUd/SztaGLiSVF4W1PrCuoNevrDnP6gnKUUmys7yQQ1Jxx1PRomaauQY6cVsqgL8Cu5lCSmxPnTGWrQVKhSGKbmtoOTpk3jVMXlNPUNUhAazYf7OToisnRsrtbEpNmxSYwGvKPpMB8av3BaFrTJ2oOUlU5nQXTy6icNYm9rb0JiY8ANtQfHvU9ItpbD40kKXp1ewsLZ5SxZO7UxM4hlDWvtdeg7rrDrNzWDMDaXa30DSem64RQutDtDd3MnlpK5axJ0eXbGrqZUlpEW+8Qpy0o5939HRx/5BR6Bv1MKPbw7v4OPEoxt3wipy0oH9Xu7KmlbKoP9eX86SPZ3369di/fWnEC0yeVACNC0NA5QFPXIEGtmV5WwvZwkqbdLb0EghqvZ+TAj0021NY7xNrdrdHvkfHzRM1BPnnmQoq9ijf3tPP+wVCip+1N3QlpS9/e2860icUM+QMsnTuNnU09nDx/WrQPKqZMYF9rL0dMLWVRTP/sbOqhrWfEluHA6MvMYFAnFaYNdYfZ3dxL5axJzJlWitej2Hyoi6NnTcLjURwTHoP7WkPbP7HEG113f1sf5ROLDROZPbupIfr5pS2N/D2cLO6ZjQ3MLZ/IrMkTeHL9QSqmTOB/qkNZP597v5FT5pfz1/D32vZ+NtV3Ul5WTFPXICVFI0GCd/a1s6ell8uWzGbBjJAvuaelhyOmltI7FOrXyD7d2dRDe9/IuHxy/aHo54iEv767jSF/kCF/kN3NPVx1ylwmTxiRyNaeIZ7ecJCrT5lr2I+RVLlzppnLMJgNKl3MUil1GvAQsA04FagBbtZa98WVuwm4CWDhwoVnHDhwICODKm95HoD9d17JoltfSFn2c+cexfevOSm6zt7/ujL60ElsXZctmR0VjdlTJ9DcPcR1p8/juc2N0QyEtXddFS0/a3IJbb3DnDRvKgtnlPHC5qZonSuWzublrc3c+/FT+dgZ86PrxK4f+f71xzZEU36ee/TMaMa2S0+czW8/VzXKxgjb/t8Klnz3ZdP9FcvcaaU0pMiDPG1iMZu+dznHf+dFhvzZxY823HYZ0yeVJNg/Y1IJPn+QniG/pfoi6TGv+OnrUaGM79MI/3nlCfzXCzui+yKe9bddxrIfrDRs59yjZ1Jz4PAoYYu0/eVHanhxS9Oo8h9aeiQvbR1Z9s3Lj+OryxfT0DnAB+5aNapsidcTrffbV57IgY4+HnmnDoA7rj2JYo+Hbz35/qh1YtPCxhNJlbvmmxcxY3IJp9z+yqjfn/zyuXz0V28brhvZrkj/ffvKE5lQ7OG7z2xN2O74Pj66YlI0u16qcm6hpMjDrh9eAYRsPOHIKdGsgJ85eyF3XHtyRrafWTmdv37pAynXjU3rGqsFmaCUqtFaV5kpaybGXQQsA36ltT4d6ANuiS+ktX5Ia12lta6qqDD2fK1g5jJ2fV2nqboiqSwBmruHwuseTkiMH6GtN/kFRWRAm3mZQ6x9sWk2U11g+/yZX7+nEm0g5oUU2Qf9BwwS3EPIA7Uq2rFsb0yd4hVCnifAa7taDX/vHUze/rr9yVOCGo25+FTC28OCYLSX4r3cWOo7BthjMGZSjeFIqtz2vmHDFwpE0gSbYWdzTzTNbDpiRbsQiD+OY1O5ZhOVfK/2cPpCecKMcB8EDmqt14W/P0FIyPOO2ZiU0VWF2RtxyWKFpp4xKaCYWSGSbB+miulryCq4PDKV0MSUx5iGkpU3d/9BZ32HXjFys1oofNIKt9a6CahXSh0fXnQJobCJo2Ty8F2yg8NocaYP0hTi1KGxRmQfZCeGRvWamX+e2VTCZDaZzcmR7YyfUN6d8afchfhUpBnMzir5V+DR8IySfcCNzplkP4YPrpjdoUmSTNkxI0HIjmRilq57LSWRslA2FUENXoOGzXjuQW3P1EU35IcR7MGUcGutNwKmguZ2Yerx7riBmGwN41CJtUvd+DZMeWdpSwiZENl1yb3YdDfcrbQ1ui5LU9liCge1xmuwthlPWmud9cleocZpqGRsekmufXIys1CJcTmjg8Pc49068VCzMA7G46VpLhgJlRj/7uSDRJnOAU9uq5lxmL38KAWe8ancYxLXCredZHNzMqGuhA9Z1CE4QsqbkzrFDWcLO8Zq6CLp1YHJde3IAz8e/YixGpZ0rXBnEClJ/kirYf3mQiXxXnM0xp3ePMEh0u07O3KmR9uK+54uQZVRWcg8rBMxwjD/S/o1R2xREuMeS7hXuDOIISc7BoxnlVi3CSweLMnqGKtuQI5I13vpxk7yV4qlx8rb3WNJPnUx/br2jBZ3PgTvNGP1UHOtcNtJpjcnIfnJQcTXveRi15i6YjPxAgpzs0qMb06mE+Jk2SKtti+4D9cKd7ZpHdOVy/Q9lFbSTgoOkabv052Uk18JpV+WLFVrOpKFb8zNKjHeJsuhEoObk2N9HNsxjdKNuFa4zRA/6KzcuTc7YLPKDpjFJbmQnHQHo52zShI80gzjDclOJna9vMJwvZgVFeNzempe0zA7SGELt9lyGU4HTFVXJtMVhdyQ1uO24Sad9Scnk1Vkpi17/EbDUIkN9boFo7DPWL2iMPvkZM4xk+hmU30nr8RkbtvW2M2caaX85d26USkvjXZefFrRSPbAWDoHhtnVPHrlSGrZ92o72BmTzCaSwCnCgfY+9ibZhjU7W+nsH6bmQGISm5+v2m24jl00pUlEZZZ393dw1qIZttQVoX94dHKoZzYeMiwXm63RiPUG/RpLb1wSrL4hP89vbmSPQdrc+KETDGpe393KS1tS27B2dyv+wMjaB9r7oml2YzGTkKu1Z4gXDba5zSA1byyxicCaugZHpUUGuPPF7Zx/7Ky07f/xrVpaeuwZN04RCGqq4xKCATR1D3DAIOWsWdLdA9jZ1EN9Rz+LKialLGc3adO6ZkJVVZWurq7OaF23po60m8qZZdS29+fbDNew/84r+eeHa3jF4AQqCOn4xqXHcf/fd9le753XncytT21OW+7IqaU0hU/MbknrOmY5amZOXuRjiIj2aLQenX5XEKywy+BqyQ6MXhBhRJPB1ZSTjGvhHqvxr0JEMz6f7BOETBjfwj2mbs0UNlprebJPEEwyvoVbdNs1aGQWjiCYZVwLt+AuJJuiIJhjXAu3eNzuQWuJcQuCWca1cAvuQaNFuAXBJONauCXBjntIlSdbEITRjG/hzrcBwijkBS2CYI5xLdyZ5isR7CcU4xblFgQzjHPhzrcFQgTD93sKQr5xqUaYFm6llFcptUEp9ZyTBuUSiXG7DFFuQTCFFY/7ZmC7U4bkA/G43YPW8k5EwX24VSJMCbdSaj5wFfBbJ4352avOpjSNp6NvOKftCcnpHPCxp6U332YIBcrz7zc6Um+dS5PBmfW4fwJ8C0j6/myl1E1KqWqlVHVra2tGxty30v60jEJhcN5dq/JtgiAk8NLW1HnX80Va4VZKXQ20aK1rUpXTWj+kta7SWldVVFTYZqAgCIIwGjMe93nAR5RStcBjwHKl1COOWiUIgiAkJa1wa61v1VrP11pXAtcDq7TWNzhumSAIgmDIuJ7HLQiCUIhYelmw1noNsMYRSwRBEARTiMctCIJQYIhwC4IgFBgi3IIgCAWGCLcgCEKBIcItCIJQYIhwC4IgFBgi3IIgCAWGCLcgCEKBIcItCIJQYIhwC4IgFBgi3GOQkqL87NY7rj0pL+0KwnhDhHsMkq93aSp5aaQg5AQR7jFIvt6lqV37hj5BGFuIcI9BAvIWZEEY04hwC7aRpwiNIIw7RLgF2xDdFoTcIMItCIJQYIhwC4IgFBgi3IJ9SJBbEHKCCLcgCEKBIcItCIJQYKQVbqXUAqXUaqXUdqXUVqXUzbkwTCg8JFAiCLmhyEQZP/BvWuv1SqkpQI1SaqXWepvDtgkFhoS4BSE3pPW4tdaNWuv14c89wHZgntOGCYIgCMZYinErpSqB04F1Br/dpJSqVkpVt7a22mNdAXLesTMtr3N0xSRbbbjxvEpb6zNDWYmXU+ZPs7ROsTd3SamOsbmPxyPlZcWOt3HDOQsdb2MsoMxmklNKTQZeA+7QWj+VqmxVVZWurq7OyKAhf4ASr4fhQBCPUhR7PQz7gwS1prTYiy8QpMgTOuAHfAGKPJ5oGlOtNf6gxqMUvkCQYq8HBfQO+5lQ5Ileyns9oXojZXyBIP6AZmKJF4C+IT8Ti70Mh9vSQLHXgz8QxB/UeD0Kj1IM+QOUlYSiTZG6ItsQsd0XCBIIaoo8ioDWTCjyMuQPRG0pLfZGtz3WnqAObQdAUGuKPB78wSATirzROkuLQ3VNKPLiDwQZ9AcpLfJQ5E08Hw/6AigVymNSWuQlqDVBDUoRtTvSvj8QBIjWE+mfYq/CF9AUeRVBrfEHNGUlXpQaEeDYbY+sG9Qar1KjtifSp5Eyw/4gpcVevJ6RunxhOwLhPo8Q6SNFKKFWsVcRCOrovinyKJRSBIMaXzCI1kT7qtjjwR/eH55wnYO+QLT+2H0WCIbsnFTiRWsY9AfwBzUTi71R2+NtLfZ66BvyU+RV0X7QWtM75MejFBOLvdF2g0HNcCAYHQORsa3C4xegyKPoHw4waUIR/vAxEdl3HhUay8OBIMUeT2ichO2K1K1UKGujRxEes0E8YdMnFHkZ9AWYUORBhcezVym8MTYUez0M+kLHpMejDMd5MDyuAQaGAxR5VUL/xLYDoWM1EBzZjsh40FrjC4QOjtj0xJFxHuknXyBIiTc01iPHU5FHRY//iF7Ej8dYuyP7a9AXoMgT2u7+sP2R42zYH+rzIu9Iea01Q+HxGnvMdg/4mFDsZfIEMxHoRJRSNVrrKlNlzQi3UqoYeA54WWt9X7ry2Qi3IAjCeMSKcJuZVaKA3wHbzYi2IAiC4CxmYtznAZ8FliulNob/rnTYLkEQBCEJpmPclipVqhU4kOHqs4A2G82xC7HLGmKXNcQua4xFu47SWleYKeiIcGeDUqrabJwnl4hd1hC7rCF2WWO82yWPvAuCIBQYItyCIAgFhhuF+6F8G5AEscsaYpc1xC5rjGu7XBfjFgRBEFLjRo9bEARBSIFrhFsp9SGl1E6l1B6l1C05btswda1S6nal1CGj+etKqVvDtu5USq1w0LZapdTmcPvV4WUzlFIrlVK7w/+nh5crpdTPwna9r5Ra5pBNx8f0yUalVLdS6uv56i+l1H8rpVqUUltillnuI6XU58LldyulPueQXT9SSu0It/20Uqo8vLxSKTUQ03cPxqxzRngM7AnbnlWSlyR2Wd53dh+zSex6PMamWqXUxvDynPRXCm3I7/jSWuf9D/ACe4GjgRJgE7Akh+3PAZaFP08BdgFLgNuBbxqUXxK2cQKwKGy71yHbaoFZccvuAW4Jf74FuDv8+UrgRUAB5wDrcrTvmoCj8tVfwIXAMmBLpn0EzAD2hf9PD3+e7oBdlwNF4c93x9hVGVsurp53gXPDNr8IXOGAXZb2nRPHrJFdcb//GPhuLvsrhTbkdXy5xeM+C9ijtd6ntR4GHgOuyVXj2nrq2muAx7TWQ1rr/cAeQtuQK64B/hj+/EfgH2KW/0mHeAcoV0rNcdiWS4C9WutUD1w52l9a67VAh0GbVvpoBbBSa92htT4MrAQ+ZLddWutXtNb+8Nd3gPmp6gjbNlVr/bYOKcCfYrbFNrtSkGzf2X7MprIr7DV/AvhLqjrs7q8U2pDX8eUW4Z4H1Md8P0iecn6rxNS1Xw1f8vx35HKI3NqrgVeUUjVKqZvCy2ZrrRshNLCAI/JgV4TrGX0w5bu/Iljto3zY+E+EvLMIi5RSG5RSrymlLggvmxe2JRd2Wdl3ue6vC4BmrfXumGU57a84bcjr+HKLcBvFoHI+3UWFUtc+CXxda90N/Ao4BjgNaCR0qQa5tfc8rfUy4ArgK0qpC1OUzWk/KqVKgI8Afw0vckN/pSOZLbnuu28TervUo+FFjcBCrfXpwP8F/qyUmppDu6zuu1zv008x2kHIaX8ZaEPSoknat9Uutwj3QWBBzPf5QEMuDVCh1LVPAo/qcL5xrXWz1jqgtQ4Cv2Hk8j5n9mqtG8L/W4CnwzY0R0Ig4f8tubYrzBXAeq11c9jGvPdXDFb7KGc2hm9MXQ18Jnw5TzgU0R7+XEMofnxc2K7YcIojdmWw73LZX0XAdcDjMfbmrL+MtIE8jy+3CPd7wGKl1KKwF3c98GyuGg/HzxJS18bFh68FIne7nwWuV0pNUEotAhYTuiFit12TVOg9nyilJhG6sbUl3H7krvTngGdi7PrH8J3tc4CuyOWcQ4zygvLdX3FY7aOXgcuVUtPDYYLLw8tsRSn1IeA/gI9orftjllcopbzhz0cT6qN9Ydt6lFLnhMfpP8Zsi512Wd13uTxmLwV2aK2jIZBc9VcybSDf4yvTu5p2/xG6G7uL0Jnz2zlu+3xCly3vAxvDf1cCDwObw8ufBebErPPtsK07yfIufwq7jiZ0t34TsDXSL8BM4FVgd/j/jPByBfwybNdmoMrBPisD2oFpMcvy0l+ETh6NgI+QZ/OFTPqIUMx5T/jvRofs2kMo1hkZZw+Gy340vI83AeuBD8fUU0VISPcCvyD84JzNdlned3Yfs0Z2hZf/AfhSXNmc9BfJtSGv40uenBQEQSgw3BIqEQRBEEwiwi0IglBgiHALgiAUGCLcgiAIBYYItyAIQoEhwi0IglBgiHALgiAUGCLcgiAIBcb/Bykqwdz/xyoeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sizes1, sizes2 = trace.get_sampler_stats('depth', combine=False)\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, sharey=True)\n",
    "ax1.plot(sizes1)\n",
    "ax2.plot(sizes2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1c21520438>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADvBJREFUeJzt3X+s3XV9x/HnGyq4H0iVFtO0xctm3ejMJqwBEpPpwCyAjuIEA+pWTLcGg84Fl8nmkjm3ZbglombEpAKxkPFrbAkdwS2utCGate4iBQSCFMagK6HXUXCG4MS998f5VO5ub3u/9/b86pvnI7m53+/n++k5r3vOva/7vd/z/Z5GZiJJquuoUQeQJA2WRS9JxVn0klScRS9JxVn0klScRS9JxVn0klScRS9JxVn0klTcolEHAFiyZElOTEyMOoYkHVHuvffe72bm0rnmjUXRT0xMMDk5OeoYknREiYj/6DLPQzeSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVNxYXBkrSUeym3Y8teB/+4EzTupjktm5Ry9JxVn0klScRS9JxVn0klScRS9JxVn0klScRS9JxVn0klScRS9JxVn0klScRS9JxVn0klScRS9JxVn0klScRS9JxVn0klScRS9JxVn0klRc56KPiKMj4r6IuLOtnxwROyLisYi4NSKOaePHtvVdbfvEYKJLkrqYzx79x4FHpq1/Frg6M1cB+4D1bXw9sC8z3wxc3eZJkkakU9FHxArg3cC1bT2As4Db25RNwAVteW1bp20/u82XJI1A1z36zwN/APxvWz8BeD4zX27ru4HlbXk58DRA2/5Cmy9JGoE5iz4i3gPszcx7pw/PMjU7bJt+uxsiYjIiJqempjqFlSTNX5c9+rcD50fEk8At9A7ZfB5YHBGL2pwVwJ62vBtYCdC2Hw88N/NGM3NjZq7JzDVLly49rC9CknRwcxZ9Zv5hZq7IzAngYuDuzPwgsBW4sE1bB9zRlje3ddr2uzPzgD16SdJwHM559J8EroiIXfSOwV/Xxq8DTmjjVwBXHl5ESdLhWDT3lFdk5jZgW1t+Ajh9ljkvARf1IZskqQ+8MlaSirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16Sipuz6CPitRHxzYi4PyIeiog/beMnR8SOiHgsIm6NiGPa+LFtfVfbPjHYL0GSdChd9uh/AJyVmb8EvA04JyLOBD4LXJ2Zq4B9wPo2fz2wLzPfDFzd5kmSRmTOos+e77fV17SPBM4Cbm/jm4AL2vLatk7bfnZERN8SS5LmpdMx+og4OiJ2AnuBrwGPA89n5sttym5geVteDjwN0La/AJzQz9CSpO46FX1m/igz3wasAE4HTpltWvs82957zhyIiA0RMRkRk1NTU13zSpLmaV5n3WTm88A24ExgcUQsaptWAHva8m5gJUDbfjzw3Cy3tTEz12TmmqVLly4svSRpTl3OulkaEYvb8k8A7wIeAbYCF7Zp64A72vLmtk7bfndmHrBHL0kajkVzT2EZsCkijqb3i+G2zLwzIh4GbomIPwfuA65r868DboyIXfT25C8eQG5JUkdzFn1mPgCcOsv4E/SO188cfwm4qC/pJEmHzStjJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJak4i16SirPoJam4RaMOIEnj4qYdT406wkC4Ry9JxVn0klScRS9Jxc1Z9BGxMiK2RsQjEfFQRHy8jb8hIr4WEY+1z69v4xERX4yIXRHxQEScNugvQpJ0cF326F8GPpGZpwBnApdHxGrgSmBLZq4CtrR1gHOBVe1jA/ClvqeWJHU2Z9Fn5jOZ+a22/N/AI8ByYC2wqU3bBFzQltcCN2TPdmBxRCzre3JJUifzOkYfERPAqcAO4I2Z+Qz0fhkAJ7Zpy4Gnp/2z3W1MkjQCnYs+In4a+Hvg9zLze4eaOstYznJ7GyJiMiImp6amusaQJM1TpwumIuI19Er+bzPzH9rwsxGxLDOfaYdm9rbx3cDKaf98BbBn5m1m5kZgI8CaNWsO+EUgSQtR9aKnw9HlrJsArgMeyczPTdu0GVjXltcBd0wb/6129s2ZwAv7D/FIkoavyx7924HfBB6MiJ1t7I+Aq4DbImI98BRwUdt2F3AesAt4EfhwXxNLkuZlzqLPzK8z+3F3gLNnmZ/A5YeZS5LUJ14ZK0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVJxFL0nFWfSSVFyX/zNWkobuph1PjTpCGRa9pIGxrMeDh24kqTiLXpKKs+glqTiLXpKKs+glqTiLXpKKs+glqTjPo5c0J8+HP7K5Ry9JxVn0klScRS9JxVn0klScRS9JxXnWjfQq4Zkzr17u0UtScRa9JBVn0UtScXMeo4+I64H3AHsz861t7A3ArcAE8CTw/szcFxEBfAE4D3gRuDQzvzWY6NKrk8faNV9d9ui/ApwzY+xKYEtmrgK2tHWAc4FV7WMD8KX+xJQkLdScRZ+Z9wDPzRheC2xqy5uAC6aN35A924HFEbGsX2ElSfO30GP0b8zMZwDa5xPb+HLg6WnzdrcxSdKI9PvF2JhlLGedGLEhIiYjYnJqaqrPMSRJ+y206J/df0imfd7bxncDK6fNWwHsme0GMnNjZq7JzDVLly5dYAxJ0lwWemXsZmAdcFX7fMe08Y9GxC3AGcAL+w/xSHqFZ85omLqcXnkz8E5gSUTsBv6EXsHfFhHrgaeAi9r0u+idWrmL3umVHx5AZknSPMxZ9Jl5yUE2nT3L3AQuP9xQkqT+8cpYSSrOopek4nybYukw+KKqjgTu0UtScRa9JBVn0UtScRa9JBVn0UtScZ51o1c9z5xRde7RS1JxFr0kFWfRS1JxFr0kFWfRS1JxFr0kFWfRS1JxFr0kFecFUyrDC5+k2blHL0nFuUevseJeudR/7tFLUnEWvSQVZ9FLUnEWvSQVZ9FLUnEWvSQV5+mVGghPk5TGh0Wvg7KspRo8dCNJxVn0klScRS9JxXmM/gjh8XJJC+UevSQVZ9FLUnEWvSQVZ9FLUnGv6hdjF/oC5wfOOGno9ylJCzWQPfqIOCciHo2IXRFx5SDuQ5LUTd+LPiKOBq4BzgVWA5dExOp+348kqZtBHLo5HdiVmU8ARMQtwFrg4QHc10h4+EXSkWQQRb8ceHra+m7gjAHcD2DpStJcBlH0MctYHjApYgOwoa1+PyIenTFlCfDdPmfrl3HOBuOdz2wLN875zLZAHzy8fG/qMmkQRb8bWDltfQWwZ+akzNwIbDzYjUTEZGau6X+8wzfO2WC885lt4cY5n9kWbhj5BnHWzb8BqyLi5Ig4BrgY2DyA+5EkddD3PfrMfDkiPgr8M3A0cH1mPtTv+5EkdTOQC6Yy8y7grsO8mYMe1hkD45wNxjuf2RZunPOZbeEGni8yD3idVJJUiO91I0nFjbzo53q7hIi4IiIejogHImJLRHQ6nWhI2S6LiAcjYmdEfH2YVwB3fZuJiLgwIjIihnrWQYfH7tKImGqP3c6I+O1xydbmvL993z0UETeNS7aIuHraY/adiHh+WNk65jspIrZGxH3tZ/a8Mcr2ptYhD0TEtohYMcRs10fE3oj49kG2R0R8sWV/ICJO62uAzBzZB70Xax8HfgY4BrgfWD1jzq8CP9mWPwLcOkbZXjdt+Xzgn8YlW5t3HHAPsB1YM2bP66XA34zp99wq4D7g9W39xHHJNmP+x+id7DBOj91G4CNteTXw5Bhl+ztgXVs+C7hxiI/drwCnAd8+yPbzgK/Suw7pTGBHP+9/1Hv0P367hMz8H2D/2yX8WGZuzcwX2+p2euflj0u2701b/SlmuTBsVNmaPwP+CnhpSLn265pvFLpk+x3gmszcB5CZe8co23SXADcPJVlPl3wJvK4tH88s19CMMNtqYEtb3jrL9oHJzHuA5w4xZS1wQ/ZsBxZHxLJ+3f+oi362t0tYfoj56+n91huGTtki4vKIeJxeof7uuGSLiFOBlZl555AyTdf1eX1f+zP19ohYOcv2QeiS7S3AWyLiGxGxPSLOGaNsQO8wBHAycPcQcu3XJd+ngQ9FxG56Z959bDjROmW7H3hfW34vcFxEnDCEbF3MtwvnZdRF3+ntEgAi4kPAGuCvB5po2l3OMnZAtsy8JjN/Fvgk8McDT9VzyGwRcRRwNfCJIeWZqctj94/ARGb+IvAvwKaBp+rpkm0RvcM376S313xtRCwecC6Yx88DvQsRb8/MHw0wz0xd8l0CfCUzV9A7HHFj+34ctC7Zfh94R0TcB7wD+E/g5UEH62g+z/28jbroO71dQkS8C/gUcH5m/mCcsk1zC3DBQBO9Yq5sxwFvBbZFxJP0jvltHuILsnM+dpn5X9Oeyy8Dvzwu2dqcOzLzh5n578Cj9Ip/HLLtdzHDPWwD3fKtB24DyMx/BV5L771cRp4tM/dk5m9k5qn0+oTMfGEI2bqYb9/Mz7BejDjICxCLgCfo/Qm6/wWUX5gx51R6L7KsGsNsq6Yt/zowOS7ZZszfxnBfjO3y2C2btvxeYPsYZTsH2NSWl9D7k/qEccjW5v0c8CTtOpgxe16/Clzalk+hV1YDz9kx2xLgqLb8F8Bnhvz4TXDwF2Pfzf9/Mfabfb3vYX6hB/kCzwO+08r8U23sM/T23qH3Z/2zwM72sXmMsn0BeKjl2nqosh12thlzh1r0HR+7v2yP3f3tsfv5McoWwOfo/R8KDwIXj0u2tv5p4KphPp/zeOxWA99oz+tO4NfGKNuFwGNtzrXAsUPMdjPwDPBDenvv64HLgMumfc9d07I/2O+fV6+MlaTiRn2MXpI0YBa9JBVn0UtScRa9JBVn0UtScRa9JBVn0UtScRa9JBX3f5nWMf7LR+lPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "accept = trace.get_sampler_stats('mean_tree_accept', burn=1000)\n",
    "sb.distplot(accept, kde=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8156558694178427"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "accept.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the index of all diverging transitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([], dtype=int64),)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace['diverging'].nonzero()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is often useful to compare the overall distribution of the\n",
    "energy levels with the change of energy between successive samples.\n",
    "Ideally, they should be very similar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c21831e10>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4XNWZ+PHvO6Mp6r25yN24F9zAgEMoDhBKWDrZUDZZQhJ2WZLdLPllQwibvklIgxAglEAINcUB00MgMca4Ytwty7Ikq/c60pTz++OObVmojGVJd6R5P8+jRzP3njvzXtl6defcc94jxhiUUkrFBofdASillBo5mvSVUiqGaNJXSqkYoklfKaViiCZ9pZSKIZr0lVIqhmjSV0qpGKJJXymlYogmfaWUiiFxdgfQU1ZWlpk8ebLdYSil1KiyefPmWmNM9kDtoi7pT548mU2bNtkdhlJKjSoiciiSdtq9o5RSMUSTvlJKxRBN+kopFUMi6tMXkQuAnwFO4GFjzPd77F8F/BRYAFxrjHm+x/4UYDfwR2PMbUMRuFIq+vj9fsrKyvD5fHaHMmZ5vV4mTJiAy+Ua1PEDJn0RcQL3AecDZcBGEVljjNnVrVkJcBPwn328zP8Cbw8qQqXUqFFWVkZycjKTJ09GROwOZ8wxxlBXV0dZWRlTpkwZ1GtE0r2zHCg0xhQZY7qAp4HLegRSbIzZDoR6HiwiS4Bc4LVBRaiUGjV8Ph+ZmZma8IeJiJCZmXlSn6QiSfrjgdJuz8vC2wYkIg7gx8B/DdDuFhHZJCKbampqInlppVSU0oQ/vE725xtJ0u/tHSJdY/GLwFpjTGl/jYwxDxpjlhpjlmZnDzi3QCml1CBFciO3DJjY7fkEoDzC1z8dOEtEvggkAW4RaTXG3HliYSqlRqOnNpQM6etdv6JgSF8vFkWS9DcCM0RkCnAYuBa4PpIXN8Z8+shjEbkJWKoJX9mlrwSkiUT1JxAIEBcXdcULBm3A7h1jTAC4DXgVa9jls8aYnSJyj4hcCiAiy0SkDLgK+LWI7BzOoJVSqi9PPvkky5cvZ9GiRXz+858nGAySlJTE17/+dRYuXMhpp51GVVUVADU1NVxxxRUsW7aMZcuWsW7dOgDuvvtubrnlFlavXs0NN9xAe3s7V199NQsWLOCaa65hxYoVbNq0id/85jfccccdR9/7oYce4stf/rIt5x2piCZnGWPWGmNmGmOmGWO+E952lzFmTfjxRmPMBGNMojEm0xgzt5fXeEzH6CulhtPu3bt55plnWLduHdu2bcPpdPK73/2OtrY2TjvtND744ANWrVrFQw89BMDtt9/OHXfcwcaNG3nhhRf43Oc+d/S1Nm/ezJ///Geeeuop7r//ftLT09m+fTvf+MY32Lx5MwDXXnsta9aswe/3A/Doo49y8803j/yJn4Cx85lFKRXz3nzzTTZv3syyZcsA6OjoICcnB7fbzcUXXwzAkiVLeP311wF444032LXr2JSj5uZmWlpaALj00kuJj48H4B//+Ae33347APPmzWPBggUAJCYmcs455/Diiy8ye/Zs/H4/8+fPH5mTHSRN+kqpMcMYw4033sj3vve947b/6Ec/OjrU0el0EggEAAiFQqxfv/5ocu8uMTHxuNfty+c+9zm++93vMmvWrKi/ygetvaOUGkPOPfdcnn/+eaqrqwGor6/n0KG+Kw6vXr2aX/7yl0efb9u2rdd2Z555Js8++ywAu3bt4sMPPzy6b8WKFZSWlvLUU09x3XXXDcVpDCu90ldKDZuRHhk1Z84cvv3tb7N69WpCoRAul4v77ruvz/Y///nP+dKXvsSCBQsIBAKsWrWKBx544CPtvvjFL3LjjTeyYMECFi9ezIIFC0hNTT26/+qrr2bbtm2kp6cPy3kNJenvY4sdli5danQRFTUcdMjm8Nu9ezezZ8+2O4whFwwG8fv9eL1eDhw4wLnnnsu+fftwu90AXHzxxdxxxx2ce+65IxJPbz9nEdlsjFk60LF6pa/GpKGeFKRiW3t7Ox//+Mfx+/0YY/jVr36F2+2msbGR5cuXs3DhwhFL+CdLk75SSg0gOTm512Vc09LS2Ldvnw0RDZ7eyFVKqRiiSV8ppWKIJn2llIohmvSVUiqG6I1cpXrR2+gfHdo5CJseHdrXWxr9M14j8dhjj7Fp0yZ++ctf8sADD5CQkMANN9zAnj17uPbaaxERnn/+eV566SV+9atfceqpp/K73/1uSN5bk76KCdNKnvvItgMFV9kQiRqNhrO88q233nr08Z/+9Ccuu+wyvvWtbwFw//338/LLLw96PdzeaNJXSo0pTz75JD//+c/p6upixYoV3H///TidTpKSkrj99tt58cUXiY+P589//jO5ubnU1NRw6623UlJifbr76U9/yhlnnMHdd99NeXk5xcXFZGVl8fDDD3PTTTexZ88eZs+eTXFxMffddx8ffPABO3bs4N577wWs8sq7d+/mJz/5yXFxPfroo3zve98jPz+fmTNn4vF4AKuMc1JSEnPmzOGnP/0pTqeTd955h1NOOYWioiIuvfRS/uVf/uW4Es4nQ/v0lVJjRl+llQFbyytXVFTwzW9+k3Xr1vH6668fV9nziIsuuohbb72VO+64g7feeosHHniAcePG8dZbbw1Zwge90ldjnDGGLSWNPL63AF/IgccR4spxtSxMabc7NDUM+iqtDNhaXnnDhg2cffbZHFkD/JprrrFtUpcmfTVmNfv8vLC5jP3VrRTEO8hx+yn1efju/gJOT2/mE+NCuON6fNgN33icVlJ/dJP2/Y8efZVWBnC5XLaWVz7y3nbT7h01JrV1BvjNPw5SXNfGJQvy+cHsYv5r+mF+NOcgV+XX8F5DMk9uOEQgGLI7VDWETrS0MoxMeeUVK1bwt7/9jbq6Ovx+P88999GBBSNFr/TVmNPaGeDx9cU0tHVx0xmTmZqVhCM8AtPtMFw5ro4sd4BfHRKe3ljK9SsKiHPq9c+wGOEhln2VVp40aVKfx4xEeeX8/HzuvvtuTj/9dPLz8zn11FMJBoNDc9InSEsrqzElGDJ87vGNvL2vhk+vmMTs/BSg9yGbT/g/zovbK7hp5WTuvjS8rHO4e2fDwY927+g4/YGN1dLKEF3llU+mtHJElzcicoGI7BWRQhG5s5f9q0Rki4gEROTKbtsXich6EdkpIttF5JpI3k+pwfrhK3t4a28NlywcdzTh92XltCzOmJbJY+8W8+R7/XcBKNXe3s6ZZ57JwoULufzyy48rrzxz5kzi4+NHRXnlAbt3RMQJ3AecD5QBG0VkjTGm+5ijEuAm4D97HN4O3GCM2S8i44DNIvKqMaZxSKJXqpvnN5fx63eK+Mxpk/pP+MaQ3rKH5R9+k4sbNuOJr6TzZSedf0vCk5gKyXmkB/NoSppOyOE69inBmXH864yR2aEqMmOlvHIkffrLgUJjTBGAiDwNXAYcTfrGmOLwvuPuihlj9nV7XC4i1UA2oElfDam/76/hzhe2c8b0TO66ZA7PbSrrtZ3L38zU8r+Q1nqArrhkqjOWkDPvYt7aXoa/vYlPOA+R3Pg+M4OdBCWOhpTZlOacQ5c7tdfXUx9ljImakSpj0cl2yUeS9McDpd2elwErTvSNRGQ54AYOnOixSvVnV3kzX3hyC9NzkvjVPy/B1cdN2aS2Ek4peRoxAYrzLmT9ou9hxMn1Kwo462NdXPPr9dxV18yvT68nq2YjGc27yW7cRnrzbsqzzoRJl4DDOcJnN7p4vV7q6urIzMzUxD8MjDHU1dXh9XoH/RqRJP3e/uVO6E+NiOQDTwA3GmM+MkZORG4BbgEoKNCbZcoSyZKHje1dPL6+mGRvHI/evIwUr6vXdu6uRmaWPos/LpG9BdfS6cnEyLEEnpHo5nefW8EN973KTesy+cyEU7kgfwrlWSspqHqDiTV/gy31cOoN4NBBb32ZMGECZWVl1NTU2B3KmOX1epkwYcKgj4/kf28ZMLHb8wlAeaRvICIpwEvA/xhj3uutjTHmQeBBsEbvRPraKrZ1dAV57N1iOvxBnr91JfmpH51cA+AMtDOz9BnEBNkXTvi9yUnx8vzZDdyxMYXHSnPZ2xrP5wocdE28kpa695lc+QoN7/ya/ROvwjjiOBAs0RE9PbhcriEtDqaGXiSjdzYCM0Rkioi4gWuBNZG8eLj9H4HfGmPsm42gxpxgyPC79w9R19rFrz+zhFPykvtsu2zXd0jwVVM44Qp8fST8I5Jchl+f3sS146p5vyGZr+6aQmGbl6rM5RzMv4j01v1MPfxniLKhzkpFasCkb4wJALcBrwK7gWeNMTtF5B4RuRRARJaJSBlwFfBrEdkZPvxqYBVwk4hsC38tGpYzUTFl7YcVFNW08anF41k5LavPdlkN25h6eA3lWWfQlDw9otd2CFyeX889sw7hFMO39hawsTGJ6oyllOacTVbzTjKbdgzVqSg1oiLqnDTGrAXW9th2V7fHG7G6fXoe9yTw5EnGqNRxNh+qZ31RHWdOz2LJpI/OfjzKGE7d/X+0e7KtG7E99DkUM2x6oo9vzzrEDw9M4McHxnPblArOzDqTtNYDTK5Yy44ZXwC0e0eNLjr3XI0qTR1+XtxewdSsRD4xN6/ftpMqXiaraTvbZ/4bIad7UO+X6gpy18wSZiV18OviPEp8Xg6M/xSC4bTt/6PdPGrU0aSvRpWXPqwgGDJcvng8TkffQwIdwU4W7f0p9SmzKRp/2Um9p8dh+I+ph0lwhvjJgfE0ODMpyT2fvPr3Ye/agV9AqSiiSV+NGnsrW9hxuImzT8khM8nTb9vJ5WtJ9FWwbeZ/gJz8f/M0V5Dbpx6mqtPNM4ezqE5fTFPiZHjzHgjZUzhLqcHQAcdqVDDG8MrOCrKS3KyaYd247bVPfunNYAyzin9LQ/IMKrNOH7IY5iR3cHZWE2/UpnFJXj3bZ/47Z239Mmx/BhZdP2Tvo9Rw0it9NSoU1rRS1dzJx2bmDFwGuegt0loL2Tv5MzDEs0KvyK8F4PnyLEpzz4Nxi+Gt70Kgc0jfR6nhoklfjQrrCmtJ8sSxcEIENXDW30+HO5Pi/E8OeRxZ7gDnZzfydl0qtW1dcO5d0FQK254a8vdSajho0ldRr7rZx76qVlZMzRj4Kr9mLxS+zv6CawY9Ymcgn8qrI04MH2zdAPXFkDIB3v4+bHxkWN5PqaGkSV9FvXeL6ohzCCum9D+bFrAWQXG62V8wfEs3pLmCnJbewrv1KXQEBaasgpZKqNs/bO+p1FDRpK+iWiAUYntZI/PHp5Lk6X3cwYaD9Ww4WM/7RTX4tj5DSfbH6PT0PuFqqHw8q5GOkJO1h71Wv747CQ6+PazvqdRQ0NE7KqodqG7D5w8xf/zAffmpLYV4u+rpcGf2ujziUJqd1EGep4tnD3q5YpIPJq2E/a9DfRFkTB3W91bqZOiVvopqO8ub8MQ5mJ6TNGDb7MYP8DsTaEqaNuxxicDZmU1sqHVT3OqESWdYG99/eNjfW6mToUlfRa1gyLCroplZeckD3sB1BjpIa91Hber84+rkD6ePZTbhwPBCsRe8qZA33xqzH/SPyPsrNRjavaOi1sHaNtq7gsyLoGsns3kHDhOiNm3hCERmyXAHWJHt55VyD1+Z1wYTlkHFB/DK1yB37rGGupauiiJ6pa+i1o7yJlxOYUZO37Xyj8hq3E6bJ5f2+P6LsA211eM62d8cR1GLE7JngzsRyjaOaAxKnQhN+ioqGWPYU9HMzNxk3HH9/zd1+ZtJ7jhMfeqcEYrumPPHWTNxXy/3WOvnjjsVqnaAv33EY1EqEpr0VVQ6UNNGsy8Q0VV+RvNeAOpTZg93WB8xITHE3DQ/r5WHJ4KNXwqhgNXNo1QU0qSvotL6A1aNm2nZiQO2TW/ZQ4c7C5+n7xW0htPqcZ1sqXNR4xNIK4DEbCjbZEssSg1Ek76KSusK60iLd5GR2H8phbhAOyltxdSnnDJCkX3U6nGdGIQ3KzzWsM0Jy6D+AHQ02haTUn3R0Tsq6gRDhvVFdUzPTkIGqJKZ1rofwdBwEl07Gw7WD/pYgFmpQSYkBHmj3MO1U3yQt8BaXKVqB0z+6DKNStlJr/RV1Nld0UxTh5+pkXTtNO+h05VCmzd/BCLrnQisyutifY0LfwhIyrW6eCo/tC0mpfqiSV9FnXWFR/rz+5+F6wh1kdZ6gIbkWUNeN/9ErcrtpC3gYGudy4olb4FVgE1H8agoE1HSF5ELRGSviBSKyJ297F8lIltEJCAiV/bYd6OI7A9/3ThUgauxa92BOqbnJJES7+q3XUrbQRwmQEPyzBGKrG+nZ/txiuGdqvA9iLz5YEJQvdvewJTqYcCkLyJO4D7gQmAOcJ2I9BwQXQLcBDzV49gM4JvACmA58E0RST/5sNVYFQwZNhfXc9rUgatkprUcICguWhIKRiCy3h2p8LnncB3TEjp4uST8K5VWAJ4UqNxuW2xK9SaSG7nLgUJjTBGAiDwNXAbsOtLAGFMc3hfqcewngNeNMfXh/a8DFwC/P+nI1Zi0r6qFtq4gSyal09HV87/T8VLbDtCcOBnjiI7xCAtT2ni+Iov6zlYyPA7InQeHN8GGB8HZ41OLlmZQNomke2c8UNrteVl4WyRO5lgVg7aUNABwakH/HwiT2krwdjWMSEXNSC1IacMgrKvu1sUT7ILaffYGplQ3kST93u6QmQhfP6JjReQWEdkkIptqamoifGk1Fm051EhmopuCjIR+2+XXrgOgMWn6SIQVkemJPhKdwWP9+pkzwOnWfn0VVSJJ+mXAxG7PJwDlEb5+RMcaYx40xiw1xizNzs6O8KXVWLS1pIHFBekDjs/Pr30Xnyt92FfIOhEOgbnJ7bxb7cYYwBkHmdOhZo/doSl1VCRJfyMwQ0SmiIgbuBZYE+HrvwqsFpH08A3c1eFtSn1EQ1sXRbVtnDoprd92jpCf3LoNUdW1c8S85DYOtzspaQvX9M+eBe210FZrb2BKhQ2Y9I0xAeA2rGS9G3jWGLNTRO4RkUsBRGSZiJQBVwG/FpGd4WPrgf/F+sOxEbjnyE1dpXraVmqVLVg8sf/+/KyGrbiCHTRGY9JPscblr6sO37jNnmV916t9FSUiGvZgjFkLrO2x7a5ujzdidd30duwjwCMnEaOKEVtKGnA6hIUT+180Jb92PSGJozlx8sgEdgLGebrI9QZZV+3m+qk+a2ZufIaV9LUkg4oCOiNXRY0tJQ3Myksmwd3/tUhO/UbqUucRcnpGKLLIicAZOV28V+MmZMIbsmdZs3NDQbvDU0qTvooOoZDhg9ImFhf035/vDLST2bST6oylIxTZiVuZ46eu08Hepm79+oFOaCi2NS6lQKtsqijw1IYSalo6ae0M0NEV5KkNJX22zW7chsMEqM5YSkJHpIPIRtbKnC4A1lW7mZ3WAVkzQBxWF09m9N2HULFFr/RVVDjc2AHAuLT4ftvl1G8iJE5q0hePRFiDMi4hxJSkAO8emaTliof0yXozV0UFTfoqKpQ3dhDnEHKSvf22y6nfRH3KXAJx/U/estvKnC421IZLLYM1UaupTKtuKttp0ldR4XBjB/mpXpyOvidlOYMdZDZ+SHXGkhGMbHBW5vhpCzjY3hDuQc2aARioK7I1LqU06SvbhYyhvLFjwK6drMbtOE2AqsxlIxTZ4J2ebfXrH+3iSZsMDpc1ikcpG2nSV7ara+2iMxBifD9Jf1rJc5xy8HEMQnLrIaaVPDeCEZ64DI9hTpr/WPE1ZxxkTIFaTfrKXpr0le2O3MQdn97/lX5K2yHavPkEo3B8fm/OyPazpc6F78jw/Mzp0FIOXa22xqVimyZ9ZbtIbuJKKEBSx2GaEyeNYGQnZ2VOF10hYVNtuCRD5gzre90B+4JSMU+TvrJdJDdxkzrKcJggLaMo6S/P9hMn5lgXT1qBVWpZu3iUjTTpK1uFQpHdxE1pO4QBW5dGPFGJcYZFGf5jN3MdTsiYqjdzla006StbHW7soDMQIj+1/6Sf3H6Idm8eQWf/4/ijzcocPx82xNHUFf4UkzUTWqugpdLewFTM0qSvbLW7ohmAvNS+k7kj2EVye9mo6c8/slj6hoP1pAdrCSFsONqvH17pq/gf9gWoYpomfWWrvZUtAOQm9z0iJ7PpQxwmQEvC5BGKaujMSPThltCxLp7UCRDnhYNv2xuYilma9JWt9lS2kJHoxuNy9tkmp34TBmgeRf35R7gchlnhJRQBq/Ba5nQ4+Hd7A1MxS5O+stXuymbyUvrvp8+t30i7N5dgXP/9/tFqXnI7+5rjqPaFf90yp0PDQWgstTcwFZO0tLIaUd3LJvuDIQ7WtHH2KTl9tneE/GQ1fEBt2oKRCG9YzEu2iqytr3ZxWUGndTMXoPjvsOh6GyNTsUiv9JVtqps7MfR/EzejaSdxId+ouYnbmykJPlJcoWPj9ZPzrCUUtYtH2UCTvrJNZbNVfiG/n+6dnPqNALQkjN6k7xCrANtx/fpTzoKD74Ax9ganYo4mfWWbyiYfLqeQkeTus01u/SYak6ZHff38gYxzNlLW7mTNnhY2HKxnI/Oguczq21dqBEWU9EXkAhHZKyKFInJnL/s9IvJMeP8GEZkc3u4SkcdF5EMR2S0iXxva8NVoVtnsIyfZi0N6L78gIT9ZDVujej3cSB3p19/RbP3xqspcbu04+I5dIakYNWDSFxEncB9wITAHuE5E5vRo9lmgwRgzHbgX+EF4+1WAxxgzH1gCfP7IHwSlKpt8A/Tn78IV7KAqI/rr5w9kvLeLdJefD1sSAWhOnEKHJ0v79dWIi+RKfzlQaIwpMsZ0AU8Dl/VocxnwePjx88C5IiKAARJFJA6IB7qA5iGJXI1qLT4/bV3Bfodr5tRvAhgVK2UNRATmJrezsyXB6sYXsf6YFf9d+/XViIok6Y8Hug8oLgtv67WNMSYANAGZWH8A2oAKoAT4kTGm/iRjVmNAZbMP6H/kTm7DJpoSp9LpyRypsIbVvOR2mgNxlPqs2cdVmSusOjy1+2yOTMWSSJJ+bx2uPS9N+mqzHAgC44ApwFdEZOpH3kDkFhHZJCKbampqIghJjXZVTVbSz+3jSl9CAbLrt4yKpREjNT+lDYAPj/TrZ2i/vhp5kST9MmBit+cTgPK+2oS7clKBeuB64BVjjN8YUw2sAz5yV84Y86AxZqkxZml2dvaJn4UadSqbfSR74kjy9D4/ML15D65g+5i4iXtEljtAnqeLHeF+/daECZA6UZO+GlGRJP2NwAwRmSIibuBaYE2PNmuAG8OPrwT+aowxWF0654glETgN2DM0oavRbKCbuMf688dO0geYm9zG7pZ4giFjdfRPPsuquBkK2R2aihEDJv1wH/1twKvAbuBZY8xOEblHRC4NN/sNkCkihcCXgSPDOu8DkoAdWH88HjXGbB/ic1CjTDBkqG7p7LNrByCvfgNNiZPxebJGMLLhNz+5nY6Q8+i6wEw5CzrqoXqnvYGpmBFR7R1jzFpgbY9td3V77MMantnzuNbetqvYVtfaSSBk+rzSd4T8ZNdv5uD4noPERr+54fH6B2paKchIsK70wRq6mTffxshUrNAZuWrEHR2508eVfmbjdlzBDiozTxvJsEZEiivIpHgfB2parQ1pEyF9ivbrqxGjVTbViKts9uEQyOlj4ZTZRY9gEBI6DjOt5LkRjm74zUtu59XaePzBcD/+1I/Bhy9AoAvi+i5JodRQ0Ct9NeIqm3xkJXmIc/b+3y+17SBt8fkEnaOzfv5A5qW0EQgZDtVZXT1MPx+6WqB0g72BqZigSV+NuMrmvkfuxAXaSWw/THPilBGOauTMTurAIRzr4pn6MXC4oPB1ewNTMUGTvhpRPn+QxnZ/n/35OfWbcBCiaQwn/XhniAnpCRQdSfqeZCg4Dfa/YW9gKiZo0lcjqmqAm7h5de8REictCRN73T9WTMtOpKyhg2af39ow43xr2GbTYXsDU2OeJn01oiqa+q+5k1u3gZaEAozDNZJhjbhp2UkYYENRuBTV9POt74V6ta+GlyZ9NaKqmn14XQ5S4z+a1ON9VaS37BvTXTtHFGQkEOcQ1hXWWhtyZkPKeO3XV8NOk74aUZVNPnJTvEgvC6eMq7FqyzcmzxzpsEZcnNPB5MxE1h+oszaIwPTzoOhtCPrtDU6NaZr01Ygxxlgjd/rozx9X/Q5t3nw6PLFRdG9adiJ7q1qoaem0NsxYDZ3NcOhdewNTY5omfTViDjd20BkI9dqf7wh2kVf3HodzVllXvWPctJLnWBWyxuW/+/oLsOlRmPZxiPPCnhdtjk6NZZr01YjZU9EC9D5yJ6d+E65gB+XZq0Y6LNtMTfCR7ArxbnV4Fq470eri2f2iVt1Uw0aTvhoxe6uspN9bdc3xNe8QcHjG1KIpA3EInJbtZ121+9iKibMuhpZyKN9qa2xq7NKkr0bM7opm0hNceF3O43cYw7jqt6nKXD5mSy/0ZVVuF2XtTopawz+TmZ8AccKev9gbmBqzNOmrEbOnsqXXrp2UtoMkd5TFVNfOEWfnWTdx/1YZ7uJJyLBq7O/+iy6YroaFJn01Inz+IAdr28jt5SbuhKq/AnA452MjHZbtJiaGmJocOJb0weriqSuEmr32BabGLE36akQUVrcSDJler/QnVr5Obep82uPzbYjMXhsO1jMrvpn3ql08tq6YpzaUwKxPWju1i0cNA036akTsqQyP3OlxpZ/YfpjM5l2U5p1vR1hRYVFqK37j4GBtuABbyjiYsMwaxaPUENNFVNSI2LttPW5HPMvr/4Kzwdp2oOAqJlZZtWZK8s6zMTp7zU7qwOMIUVm4jYu6qsCZAUl51pX+Oz+GVV+xO0Q1huiVvhoRe5rimJkSwNlj3lVB5WvUp8ymbYxX1eyP22GYm9zOtqbEY/du8xZY3ys/tC0uNTZp0lcjYk+zk1mpgeO2JXRUktW4nZIY7to5YnFqK1Vdbso7wzd0k7IhOR8qt9sbmBpzIkr6InKBiOwVkUIRubOX/R4ReSa8f4OITO62b4GIrBeRnSLyoYj0XnhFjVm1rZ3U+D6a9I907cRyf/4RS1Kt/vxNjUnHNubNh/oiaKu1KSo1Fg2Y9EUZCyelAAAeg0lEQVTECdwHXAjMAa4TkTk9mn0WaDDGTAfuBX4QPjYOeBK41RgzFzgb0BKCMWZv+CZuz6RfUPEKDckzaEmcbENU0SXTHWBKgo/NxyX9BYCBvWtti0uNPZFc6S8HCo0xRcaYLuBp4LIebS4DHg8/fh44V6zauauB7caYDwCMMXXGmODQhK5Giz29JH1PVwPZjR9wKP8iu8KKOktSW9nXFk9dZ/jGR8p4iM/QUTxqSEWS9McDpd2el4W39drGGBMAmoBMYCZgRORVEdkiIl89+ZDVaLOnopksT4gs77EZpplNOwA4lH+hXWFFnSVpLRiEv1Z4rA0iVhdP0VvQ2WJvcGrMiGTIZm91bnvOD++rTRxwJrAMaAfeFJHNxpg3jztY5BbgFoCCgoIIQlKjyZ7KluO7dowhs2kH1emLaUvoef0Qu6bEd5Lh8vNmhZurJlvLSpK3AA6+Da9+HcYtPtZ46c32BKlGvUiu9MuA7uPpJgDlfbUJ9+OnAvXh7W8bY2qNMe3AWuDUnm9gjHnQGLPUGLM0Ozs2FtCIFcGQYV9VC6d0S/rxndUkdNZo104PIlYXzzuVbnxHOkEzplgll3XophoikST9jcAMEZkiIm7gWmBNjzZrgBvDj68E/mqMMcCrwAIRSQj/MfgYsGtoQlejQXFdG52B0HFX+llNH2IQSvJW2xhZdFqe3kJ70MHbR2rxiANy50H1LggF+j9YqQgM2L1jjAmIyG1YCdwJPGKM2Ski9wCbjDFrgN8AT4hIIdYV/rXhYxtE5CdYfzgMsNYY89IwnYuKQkcWTpl9JOkbQ2bTThqTpjGh6s1+joxNc5LbSXeHWFvm5RPju6yNeQugdAPU7rcWUFfqJERUhsEYsxara6b7tru6PfYBV/Vx7JNYwzZVDNpb2YxDYHqKlfST2kvx+JsozTnH5siiU5zA6nGdvFTmwRcErxPImglOt9XFo0lfnSSdkauG1e7KFqZkJVrJC8hq2kFQ4mhIPsXewKLYRRM6aQ04eKcq3MXjdFnJvmoHGF1GUZ0cTfpqWO2pbGZWfor1JBQko3knDcmnEHK6+z8whq3M6SIt3MVzVO586GyGxkP2BabGBE36ati0dgYore9gVm6ytaFmD65gB3Vp8+wNLMq5HFYXzxvl3Ubx5M6xburqKB51kjTpq2FztPzCkSv9w1sIOL00JU63MarR4ZKJPloDDt48MlHLlQCZM6ykr8soqpOgSV8NmyfWW10R+ypb2HigkmDFdupS5mAczgGOVCtz/OR6g/zxULcunrz50FYDrVX2BaZGPU36athUNnfgiXOQluAivWUfTuOnLlW7diLhFPhUgY+/Vbqp9YUnvOeGf3ZablmdBE36athUNvnITfEiImQ2fkhnXAotCZPsDmvU+KdJPgJG+Etp+Go/Pg3SJmm/vjopmvTVsDDGUNnsIy/Vi7urkdTWA9SlzrVqDaiInJIaZG6anz+W9OjiaSqFpjL7AlOjmiZ9NSwqmnz4/CHyUrwUVL6GgxB1qfPtDmtU2HCw/ujXkqR6tje42N0Yvg+SF/4Z7tEa+2pwNOmrYbG7ohmAvBQvk8rX0uHJot2ba3NUo8+qzCZcEuLJonhrQ1IuJOVYi6YrNQia9NWw2FluJf3pngZyGzZTmzpPu3YGITkuxMqMZv54yEuL/8gN3flQvA46GuwNTo1KmvTVsNhV3kxmopsZNa8B6Kidk7A6u5H2oOPY8M28+WCCsO81ewNTo5ImfTUsdlU0k58Wz+Tyl6hNXUCnO8PukEat6Yk+FqT7eaIo3pqXlVYASXmwR5dRVCdOk74acs0+PyX17SyNryC9ZR/F43SxlJP1z9M62N8cx4Zal1WOYdZFUPgm+DvsDk2NMpr01ZDbHe7PP9f/N0Li5FD+BTZHNPpdMsFHqivEkwfCN3RnfRL8bVD0tr2BqVFHk74acrsqmhFCLGp8g8qs0+n0ZNod0qgXHwdXTfbxymEP1T4HTF4FnhTt4lEnTJO+GnK7yps5L7GI5M5KDo672O5wxoxPT+0gYIRnDnohzg0zzoe9L0NQl1FUkdOkr4bczvJmrveux++Mpyzn43aHM2ZMSQ5yVm4nTxXFEwiGYPal0F4Lh9bZHZoaRTTpqyHVFQhxqLqe033/oCz3XIJxCXaHNCYcmaG7IqmGig4n3/jzTpix2iq5vOtPdoenRhFN+mpI7a9u4UyzFW+wheJxn7Q7nDFnSWorue4utuzcA9ufsdbP3f4svP+w3aGpUUKTvhpSOw43cZlzHcH4LCozT7M7nDHHIXBBTgN72xL4oD4Oxi2GrlaoP2B3aGqUiCjpi8gFIrJXRApF5M5e9ntE5Jnw/g0iMrnH/gIRaRWR/xyasFW02nuojPOcW3DMvwLjiLM7nDHp41lNxDuC/GZ/grVgutMN5VvtDkuNEgMmfRFxAvcBFwJzgOtEZE6PZp8FGowx04F7gR/02H8v8PLJh6uiXdrBl3ETQBZeY3coY1a8M8Q5WU2sLfNQ0emF3LnWwio6ikdFIJIr/eVAoTGmyBjTBTwNXNajzWXA4+HHzwPniljVtUTkU0ARsHNoQlbRqjMQZFnLG9TH5UD5B0wree7olxpaF+Q0EDLw2wPxkL8IutrgoE7UUgOLJOmPB0q7PS8Lb+u1jTEmADQBmSKSCPw38K2TD1VFuwOF+1ghu6nLWqYVNYdZjsfPsrQWflvo5Z2WcQQcXorf1Ju5amCRJP3efntNhG2+BdxrjGnt9w1EbhGRTSKyqaamJoKQVDTybf49DjEkTTrV7lBiwkW5DbQFnbzTmEld6jwmVr2h5ZbVgCJJ+mXAxG7PJwDlfbURkTggFagHVgA/FJFi4D+A/ycit/V8A2PMg8aYpcaYpdnZ2Sd8EioKGMOEQ39gK7PIy9KyCyPhlMQOpiZ0sLYqnaq0RThDXbDjBbvDUlEukqS/EZghIlNExA1cC6zp0WYNcGP48ZXAX43lLGPMZGPMZOCnwHeNMb8cothVNClZT05XKZszL9aenREiAp/MbaC808O7XdNpSJ4JW5+0OywV5QZM+uE++tuAV4HdwLPGmJ0ico+IXBpu9husPvxC4MvAR4Z1qrEtsOlxWkw87TMuHbixGjKnpTWT7vLzUnUmRRMut4ZuVumYCdW3iAZSG2PWAmt7bLur22MfcNUAr3H3IOJTo4GvCdn1J/4SXMmcgjxotzug2BHngE9kN/B0eQ4bks9jieMnsOUJuPD7doemopTOyFUnb8cLOIM+ngmezaKCNLujiTnnZTfilhB/PRSAOZfCtt+Br9nusFSU0qSvTt6W31LmmkJT+nyykjx2RxNzkuNCrMpsYltpI40Lb4HOZtjyW7vDUlFKk746OWWboXwrj3V+nIwkD09tKDlaEVKNnItyGwiEDL85mA6TzoAND+gMXdUrTfrq5Gx8iJArkae7VlKQkWh3NDFrvLeLueNSeOzdYtqX3gpNpVpyWfVKk74avLY62PEHisZfQisJFGRq7Xw7nT0zhxZfgEdrZkHmdHj3F2B6zqNUsU7LIKrB2/pbCHbyl9Y5JDiDnN7wFxyNdgcVu8anx/Oxmdk8su4Q/7r6S7hfvgPWfhVye9ZHBJbePPIBqqigV/pqcEJB2PgIZEzjlebJzEzswKGTsmz3b+dMp66ti8fbV0L6ZNj7EpiQ3WGpKKJJXw3O3rXQVELbhLPY1+xkZlKH3REpYOnkDM4+JZtfvH2ItpVfhebDUPGB3WGpKKJJXw3Ou7+EtAI2uZZgEGYmatKPFl+7cDatnQF+VrUAkvOsP9ChoN1hqSihSV9FbtOj1tfr34TS92D8Et6t8eISo1f6UeSUvGSuXDKBx9aXUTvpYmirgdL37Q5LRQlN+urEFb0FrniYeBrv1bhZlOHH49BRInabVvLc0T/MX87dhpMgXy1diUmfYvXt+/UPs9Kkr05UWy1UbIeClTQbLx82xHF6jt/uqFQPefEh/nNeK3+t8vJ25jXWylr7X7M7LBUFNOmrE1P0FogDpqzi/RoXIYTTs7vsjkr14qbpHSzK8HPHvnn4xq2wllNsrbI7LGUzTfoqch2NVl/+xOXgTeXdajceh2Fxpl7pRyOnwA+XNNMWEL5QeyUBiaNh0/NsOFjPUxtK7A5P2USTvorcgTetGZ7TzwNgfY2LJZl+vE6b41J9mpka5H8WtPJW8zhe9VxAeut+Ulv22x2WspEmfRWZ5gooWQ8TlkNCJvWdwu4mF6fnaNdOtPvMtA7OSG/ijoYraIrLYlLla0hIP53FKk36KjLrfmbN7JxxPgDv1bgBOD1bk0e0E4FbJlWS4w1xZ8cNxHfVMfPQ03aHpWyiSV8NrOEQbHoEJiyDBGvR879Vukl2hViYoUl/NPA6DV+bUcr7jgWsM/OZt/9+aySWijma9NXA3rzHGrEz80LA6tZ/q9LNqtwuXPo/aNTIcgf42owyvh/4NM5gO52v3WN3SMoGWmVT9a9sM+x4Hlb9F8RbSyE+s6uNGl8Ok+Kq2HBQl+WLJgMtXjMxvosrpjl54uBqbvzgt/iW/gveiYtGKDoVDfQ6TfXNGHjt65CYDWfcfnTz1qYkABamtNkVmToJs5M7KJ3/7zSaRIqf/Df8Aa3LE0siSvoicoGI7BWRQhG5s5f9HhF5Jrx/g4hMDm8/X0Q2i8iH4e/nDG34aljteMEasXP218CTfHTz1qZEpiZ0kObSZDFaTZ80kf1z/4NZndv5/WO/wOhiKzFjwKQvIk7gPuBCYA5wnYj0XJXhs0CDMWY6cC/wg/D2WuASY8x84EbgiaEKXA2zjgZ45U4YtxiW3HR0c32nsL8tnsWpepU/2q248svUJM7gvNKf8eOXttodjhohkVzpLwcKjTFFxpgu4Gngsh5tLgMeDz9+HjhXRMQYs9UYUx7evhPwiohnKAJXw+z1u6C9Hi75GTiOzb56p8qNQTg1tdXG4NSQcDjJuvoXjJN6Et67lwffOWB3RGoERHIjdzxQ2u15GbCirzbGmICINAGZWFf6R1wBbDXGdA4+XDUiiv8BW34LK/8d8hcet+vlMg/pLj9TE3w2BaeGwrSS5zhSiWFK6kJuaXqJ1S+vIiPRw5VLJtgbnBpWkST93hbB69kB2G8bEZmL1eWzutc3ELkFuAWgoKAggpDUsPE1w5++YC21d/bxt2/aAsLfKj2cndmoSyOOIaV555LRsocfJjzJ1c/ns720kVn5KUf3X79CfyfHkki6d8qAid2eTwDK+2ojInFAKlAffj4B+CNwgzGm18+PxpgHjTFLjTFLs7OzT+wM1NBa+1/QVAaXPwjuxON2/bXCTWdIOC1dh2mOJYG4JEpzz2FpcBs3J23g9xtLOFSn92zGqkiS/kZghohMERE3cC2wpkebNVg3agGuBP5qjDEikga8BHzNGLNuqIJWw+TD52H707Dqq1DQswcP1pZ5yPYGmaWrZI051elLqU5fzH+Zx5nkaeeJ9w7R0K51lcaiAZO+MSYA3Aa8CuwGnjXG7BSRe0Tk0nCz3wCZIlIIfBk40i9wGzAd+IaIbAt/5Qz5WajBO7IE4ts/hD9/CdImWROxemjvCvBWpYcLxndq185YJML78+7GFWzn/oynCRnDk+8doisQsjsyNcQimpFrjFkLrO2x7a5uj33AVb0c923g2ycZoxpuwS4r8YsDTr0Btn50ZO1bZR58wVQuGt8J+sl/TGpOmsrO6bewYP993DNtFf+9ZxovbCnjxpWTENG/9GOFzsiNdcbAh89BSwUs/uejBdV6+sMhL7neIMu1quaYtnPqZ6lNW8Cnyn7A9TNDfHi4iV+9rUM5xxJN+rGu5D0o2wgzVkNOzzl3lqoOB29VuLlikg+nXvCNacbhYt3CHwLCHU3f59TxCfzfq3t5a0+13aGpIaJJP5Yd3gI7n4fsU2DmJ/ps9odDXkIIV03Wsflj2bSS55hW8hx5te9SnH8BWU07+E7875mdm8ztT2+lpK7d7hDVENCkH6va6+HZG8GdDIs/Y/Xn98IYeK7Yy/KsLqYka62dWNGQMpvyzNOZXfYM/5P+Gl3BENc+tJ7H3y3W9XVHOU36sSgYgBc+C62VsORmcCf12XRLXRxFrXF6lR+DSnPPozj/QlYe/AXfmbKD8kYfL22vsDssdZI06cei178BB/4KF/0I0if12/SJogSS4kJcNEGrZ8QcEd6b/20qM5Zzecl3uTv/Pd4vrmdrSYPdkamToEk/1mx+HN67H1Z8AZbc2G/T8nYHfyn1cM0UH4lxWno3FoWcbt5e8gsqss/kpoaf87/Jf+DP20rZW9lid2hqkDTpx5LidfDSV2DaObB64OkTjxXGA3DzdL2BF8uCcQm8c+rPKJx4BZ/xP89jcT/g6799jaZ2Hb47GulyibGioRie/YzVnXPlo+Ds/5++xS/8viie09KaOVxdy+GRiVJFKeOI4/2536Q+ZQ5Ld/8fD7X9O08+vJdbvvjfuOKcA7+Aihp6pR8LOlvg99eBvx3mXw07/3is/EIfniqKpyXg4OLc/tdcVTFEhMKCq3n1zOcIpU/lS/U/oOTe8zDVe+yOTJ0AvdIf6wKd8PSnoWYvLL8FkgYufdTUJdy/J4FVuZ1MTdQbuOp4OXUbyTzjJl5fv5nldX8kdP9KHNPORmashrgeayQtvdmeIFWfNOmPZUeGZh58Gz71AAQiG3Z5/54Emv3C1+a30qwDNWLatJLnet8hDs47fRk/3rKESaV/4qoDb2LKtyBzL4fc+aC1eqKWdu+MVe8/DI98Anb/BeZcHnHCL2tz8GhhAldM8jE7TSdjqb6JwFdOdbBr8g1c2XkXh7sSYdMjsPEhaK+zOzzVB036Y1GgE7Y8Doc3wSkXwdSPRXSYMXDPB8kI8OW5WkpTDUwE7lrYysULx3NO+3e533E9wdoD8M4PoXyb3eGpXmjSH2s6GuGpq6FyO8y93CqkFqE/lHh5rdzDHXPbGJegddRVZETgpukdPH5WC0+YC/lY+w8olfGw5TFriLBfZ3NHE036Y0nVLnjo49bC5guvgymRXeEDHG53cPfWJJZndfGvM3Vcvjpxp+f4eX11PedNT+K8lm/ySPAi2Pgw/ofOgzotzxwtNOmPBcZYM20fPhe62uCml2DiR5c77Et7AL64PpUQ8KOlzVo+WQ3KhoP17Cyr48LUEn4wt4S3Uy7hc/6v0FZVhO++M6l672m7Q1SAGBNd0+uXLl1qNm3aZHcYo0dtIbz4H1D8d5h8FlzxMCTn9TsGvzt/CP713VTeqXTzwMomVo87fl3UDQd1nL4avKpOF38vh+taH+dURyEvOc5m8/y7+Z/LFuHQdTeHlIhsNsYsHaidDtkcrRoOwTv/B9t+B043zL/GWsx878sRv0RHAL6yMYW/VXr43qnNH0n4Sp2sXI+fK6dATdd1vF3yDz7Z+RZTtt7ILXtv5/Qzz+GqpRNI8brsDjOm6JX+aBIKQtFbsPkxK7mLw+rGmX4+eFNO6KXK2x18fn0qOxri+PqCVua7tNCCGn7JTXuZUPk6iYFGHghczCPyT5w1dzIXLxjHmTOy8Lq0pMNgRXqlr0k/2rXWQOl7sO9VK9G311rr2C68Dk77Aux//YRezh+Cxwvj+emuRAB+tryZc8d1aTeOGjEl+au5qvZXsO13tDlTeSh4MY/4zsYXl8yyyeksnpjO3HEpFGQmkJfiJSPRrQuzR2BIk76IXAD8DHACDxtjvt9jvwf4LbAEqAOuMcYUh/d9DfgsEAT+3Rjzan/vFdNJv6sNavdBxXYo3WCtX1sfHvXgSYHMaZC3EHLnDVgwrae6TuG54niePBBPWbuTRSmt3DyxijyvVkpUI+tAwVUAZDTuYH7h/Yyv+TsBcbHNs5RXAqfy9/ZJ7A+NIxQeZ+KOc5Cb4iE9wU2yN45kj4tkbxxJ3jiSvS5SvHHsKm8myRtHRoKb9EQ3Lqd17PUrCmw7z5E2ZH36IuIE7gPOB8qAjSKyxhizq1uzzwINxpjpInIt8APgGhGZA1wLzAXGAW+IyExjTGxM9Qz6wd9hzYb1t1uPOxqt2YpHvlqroHa/leybSo8d606E9Ckw+1Lre9pEcESW6I2xkvyORhfb6uP4e5WbrXUuQggrsrr49LhyFqe06Ux5Zav6tHm8vfR+0pt2MaX8ReZUvMrSwHpwg188NLmyCSTk0iRJNAU9NPu8NLV7aOo0NAQTqAt6KQ8m0oqXFpNAPSnUm2Ta8ZDsdZGe4GZTcT0TMxKYFP7UkOiJI9HjJMEdhzvOQcgYjIEXNpdhsH53THjbJYvGkRbvIiXehXMM3XQe8EpfRE4H7jbGfCL8/GsAxpjvdWvzarjNehGJAyqBbODO7m27t+vr/Yb0St/6FwQTtPrDTRBMqJdk7INAh7Wt+/aityHkt9of+UqfbO0L+Lq177Be47jHbdZ7DcSVCPHpkJQLybmEEnMJJuUTTMgmhBA0QshAwEB7QGj1C20BBy1+oS0gtAaEWp+DKp+Dqg4nVR0OiludNPmtKx3BMDmhkyWprZyW3szEeL1Zq6KUMXi7aknqqCDBV4nL30JWXLffyUCn9X2A36sucdMsKdSTQoNJoSKQQL2x/iB04KYLl/Vl4ujChR8nXbgI4CSIg5BxEEQI4bCei4MEj5vkeA9J8R6SvB4SvW4Sw4+T4t0kJXhIjvfgdblwu124XXF43G7ru8uNx+PC44rD5XAgwrB0Vw3l6J3xQLdLUMqAnoPAj7YxxgREpAnIDG9/r8ex4yN4zxPXVgs/nW/9h+ie4IeEgNMFDpdVl97lhbh4cMVbjzubrcfeFKuN02WNqHG6jh135PHcy60++YRM/rCnnf/34gFCbUGCNRA0g/+PkOYOkezoIt3VxfLULvK9XUxK6GRqgo8Ep86uVaOACD5PNj5PNrAAgI9M6TIGMUGcoU7rK9gV/u7DFWwnLtBOXLAdV7CdRFcqGV0NzOkqwd3VgCd4EqVFDNAe/joJAWP9ITEIG8xsPhv42tE/Ag6BhRPSeObzp5/cmwwgkqTfWybq+fGgrzaRHIuI3ALcEn7aKiJ7I4iruyyg9gSPGaSTXRj6l33tOKlzODTYA4fWCP47DKuxcB56DtGhn3P4B/DJ47bsBZ69ddDv1f+C12GRJP0yYGK35xOA8j7alIW7d1KB+giPxRjzIPBgJAH3RkQ2RfKxJprpOUSPsXAeeg7RIRrPIZIyDBuBGSIyRUTcWDdm1/RoswY4ssr2lcBfjXWzYA1wrYh4RGQKMAN4f2hCV0opdaIGvNIP99HfBryKNWTzEWPMThG5B9hkjFkD/AZ4QkQKsa7wrw0fu1NEngV2AQHgSzEzckcppaJQRGMAjTFrgbU9tt3V7bEPuKqPY78DfOckYozEoLuGooieQ/QYC+eh5xAdou4com5GrlJKqeGjpZWVUiqGjNqkLyJXichOEQmJyNJu2yeLSIeIbAt/PWBnnAPp6zzC+74mIoUisldEPmFXjCdCRO4WkcPdfv4X2R1TpETkgvDPulBE7rQ7nsEQkWIR+TD8sx819UxE5BERqRaRHd22ZYjI6yKyP/w93c4YB9LHOUTd78OoTfrADuCfgHd62XfAGLMo/DX4Ua8jo9fz6FHC4gLg/nBJjNHg3m4//7UDN7dft3IjFwJzgOvC/waj0cfDP/uoGio4gMew/p93dyfwpjFmBvBm+Hk0e4yPngNE2e/DqE36xpjdxpgTncQVdfo5j8uAp40xncaYg0AhsHxko4spy4FCY0yRMaYLeBrr30CNAGPMO1gj/7q7DHg8/Phx4FMjGtQJ6uMcos6oTfoDmCIiW0XkbRE5y+5gBqm38hfDU8Ji6N0mItvDH3ej+iN5N6P5592dAV4Tkc3hme6jWa4xpgIg/D3H5ngGK6p+H6I66YvIGyKyo5ev/q7AKoACY8xi4MvAUyJyYiuMDLFBnkdEJSzsMMD5/AqYBizC+rf4sa3BRi5qf94n6AxjzKlY3VRfEpFVdgcU46Lu9yGql0s0xpw3iGM6gc7w480icgCYCdh2U2sw50GEJSzsEOn5iMhDwIvDHM5Qidqf94kwxpSHv1eLyB+xuq16u+81GlSJSL4xpkJE8oFquwM6UcaYqiOPo+X3Iaqv9AdDRLKP3PAUkalYpR+K7I1qUEZlCYvwL+cRl2PdqB4NIik3EtVEJFFEko88BlYzen7+vele3uVG4M82xjIo0fj7ENVX+v0RkcuBX2DV7X9JRLaFa/6vAu4RkQDWal23GmOi9uZKX+cxiktY/FBEFmF1jRQDn7c3nMj0VW7E5rBOVC7wR7FqtccBTxljXrE3pMiIyO+Bs4EsESkDvgl8H3hWRD4LlNDHrP9o0cc5nB1tvw86I1cppWLImOveUUop1TdN+kopFUM06SulVAzRpK+UUjFEk75SSsUQTfpKKRVDNOkrpVQM0aSvlFIx5P8D8GAHsoXJArsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "energy = trace['energy']\n",
    "energy_diff = np.diff(energy)\n",
    "sb.distplot(energy - energy.mean(), label='energy')\n",
    "sb.distplot(energy_diff, label='energy diff')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the overall distribution of energy levels has longer tails, the efficiency of the sampler will deteriorate quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple samplers\n",
    "\n",
    "If multiple samplers are used for the same model (e.g. for continuous and discrete variables), the exported values are merged or stacked along a new axis.\n",
    "\n",
    "Note that for the `model_logp` sampler statistic, only the last column (i.e. `trace.get_sampler_stat('model_logp')[-1]`) will be the overall model logp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pm.Model()\n",
    "with model:\n",
    "    mu1 = pm.Bernoulli(\"mu1\", p=0.8)\n",
    "    mu2 = pm.Normal(\"mu2\", mu=0, sigma=1, shape=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "CompoundStep\n",
      ">BinaryMetropolis: [mu1]\n",
      ">Metropolis: [mu2]\n",
      "Sampling 2 chains: 100%|██████████| 22000/22000 [00:06<00:00, 3399.17draws/s]\n",
      "/Users/twiecki/anaconda3/lib/python3.6/site-packages/mkl_fft/_numpy_fft.py:1044: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  output = mkl_fft.rfftn_numpy(a, s, axes)\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    step1 = pm.BinaryMetropolis([mu1])\n",
    "    step2 = pm.Metropolis([mu2])\n",
    "    trace = pm.sample(10000, init=None, step=[step1, step2], cores=2, tune=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'accept', 'p_jump', 'tune'}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.stat_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both samplers export `accept`, so we get one acceptance probability for each sampler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.25      , 1.5632343 ],\n",
       "       [0.25      , 0.53344289],\n",
       "       [1.        , 0.00924178],\n",
       "       ...,\n",
       "       [1.        , 0.09734532],\n",
       "       [1.        , 0.05016997],\n",
       "       [1.        , 0.41670172]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trace.get_sampler_stats('accept')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
